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Description  of  Progress 

Dexterous  Manipulation 

One  of  the  great  deficiencies  of  today’s  robots  is  their  lack  of  flexibility.  Most  industrial 
robots  are  capable  only  of  simple  and  repetitive  tasks,  such  as  spot  welding  or  spray  painting. 
There  are  two  main  reasons  for  this  deficiency.  First,  typical  end-effectors  use  a  very 
simple  two  stick  structure.  Second,  dexterous  manipulation  (manipulation  by  grippers  with 
independently  moving  fingers)  is  not  well  understood.  Effective  techniques  for  dexterous 
manipulation  would  have  application  within  a  wide  range  of  areas,  including  industrial 
assembly,  decontamination  of  nuclear  plants,  and  exploration  of  remote  environments  (e.g., 
ocean  bottoms  or  space). 

We  are  developing  a  new  strategy,  finger  tracking  [Rus90,Rus9l],  for  the  autonomous, 
manipulation  of  objects  by  multifinger  robot  hands.  Most  of  the  earlier  efforts  devoted  to 
understanding  dexterous  manipulation  have  been  devoted  to  grasping  or  to  manipulation 
for  task-specific  problems.  The  finger  tracking  paradigm  reorients  an  object  with  a  series 
of  rotations  effected  by  fine  finger  motions,  in  which  the  hand  maintains  contact  with  the 
object  at  all  times.  It  is  common  for  humans  to  reorient  an  object  using  “extra”  fingers  — 
those  not  needed  for  grasping.  Finger  tracking  captures  and  formalizes  these  ideas. 

Our  notion  of  manipulation  refers  to  the  reorientation  of  an  object  by  a  mechanical 
hand.  The  reorientation  is  accomplished  by  fine  finger  motions,  with  the  object  held  in  the 
hand  through  the  entire  process. 

Definition  1  Manipulation  is  the  reorientation  of  a  part  inside  the  grip,  while  maintaining 
the  grip. 

This  definition  implies  that  the  reorientation  is  accomplished  with  respect  to  a  system 
of  coordinates  which  is  fixed  on  the  robot  hand.  In  the  process  of  reorientation,  the  grasped 
object  undergoes  a  Euclidean  motion,  composed  of  a  translation  and  a  rotation.  We  are 
interested  in  measuring  rotations,  and  therefore  we  abstract  out  the  translational  component 
of  the  motion. 

Definition  2  Two  congruent  objects  have  equivalent  orientation  if  there  is  a  pure  transla¬ 
tion  that  takes  one  into  the  other. 

We  see  two  basic  types  of  manipulation  based  on  the  relationship  between  the  object 
and  the  hand.  In  the  first  kind,  the  object  is  passive  inside  the  grip.  The  object  is  fixed  with 
respect  to  the  fingers  involved  in  the  grip  in  the  reorientation  step.  Its  motion  follows  the 
finger  motions.  We  call  the  algorithms  which  result  from  such  an  interaction  finger  walking 
algorithms.  In  the  second  type  of  relationship,  the  object  is  active  inside  the  grip.  The 
internal  forces  of  the  grip  Eire  used  to  move  the  object  relative  to  the  grasping  fingers.  The 
manipulator  produces  motion  by  applying  forces  which  take  into  consideration  the  geometry 
of  the  object.  We  call  the  manipulation  algorithms  which  result  from  an  active  interaction 
between  the  fingers  and  the  object  finger  tracking  algorithms.  Both  finger  walking  strategies 
and  finger  tracking  strategies  are  common  in  human  manipulation.  We  have  experimented 
with  both  strategies  using  Newton,  our  simulator  for  rigid-body  dynamics.  Our  most  recent 
efforts  have  been  focused  on  finger  tracking. 
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In  order  to  make  the  idea  of  finger  tracking  precise,  let  0  be  the  object  to  be  manipulated 
and-let  H  be  the  dexterous  robot  hand.  We  assume  that  H  has  enough  fingers  for  a  good 
grasp  of  the  object.  The  hand  also  has  some  additional  fingers,  which  we  call  free  fingers. 
Some  of  the  fingers  are  used  to  constrain  the  object  by  restricting  its  d<  <jees  of  freedom, 
while  the  other  fingers  are  used  to  generate  its  motion.  In  a  typical  manipulation  problem, 
we  are  given  the  number  of  fingers  n  of  the  hand.  Fingers  from  a  subset  of  size  m  <  n  -  1 
axe  used  to  grasp  the  object,  usually  by  assigning  each  finger  to  a  separate  face.  The  size 
of  to  depends  on  the  contact  type  [MNP90,MSS87,Ngu86].  Once  the  object  is  grasped,  the 
to  fingers  stay  fixed  in  space  and  the  ob  iect  is  constrained  to  maintain  continuous  contact 
with  them.  In  addition  to  the  m  giasprig  fingers  that  constrain  the  degrees  of  freedom, 
the  object  is  also  in  contact  with  a  free  finger,  which  tracks  along  some  curve  on  a  different 
face.  This  process  causes  the  object  to  move  relative  to  the  grasping  fingers. 

The  free  finger  tracks  a  continuous  trajectory,  while  at  each  instant,  the  m  +  1  fingers 
hold  the  object  in  a  grasp  with  some  desired  property,  for  instance  equilibrium.  Using  this 
technique,  we  can  generate  the  reorientation  of  a  grasped  object  by  commanding  a  simple, 
sliding  motion  for  the  tracking  finger. 

The  most  fundamental  question  within  this  framework  is  exactly  how  to  generate  some 
desired  motion.  The  answer  should  allow  us  to  program  a  robot  to  take  an  object  from  a 
given  initial  orientation  to  a  goal  orientation,  or  to  determine  when  such  a  program  does 
not  exist. 

To  answer  this  question,  we  have  broken  the  problem  into  two  components.  The  first  is 
related  to  the  fa.ct  that  this  form  of  manipulation  is  defined  as  a  constraint  problem.  Thus, 
an  important  aspect  is  finding  an  algorithm  to  determine  the  configuration  space  for  the 
motion  of  the  object  to  be  manipulated.  The  second  component  has  to  do  with  finding 
the  manner  in  which  the  robot  must  use  its  fingers  to  generate  some  desired  trajectory  for 
the  object  to  be  manipulated.  This  involves  finding  efficient  tracking  algorithms.  We  have 
established  a  framework  in  which  to  adress  these  algorithmic  questions,  by  using  Lie  algebra 
properties. 

Some  of  our  results  are  summarized  below: 

•  The  configuration  space  for  a  polyhedral  object.  For  the  case  of  a  polyhedral  object  held 
by  a  robot  hand  with  four  frictionless  point  contacts,  we  have  obtained  an  algorithm 
to  describe  the  configuration  space  as  a  manifold  given  by  a  closed  form  equation. 
We  have  analyzed  the  properties  of  the  configuration  space,  and  have  shown  that  it  is 
diffeomorphic  to  the  rotation  group  50(3).  Furthermore,  we  have  shown  that  for  this 
configuration  space,  the  vertices  of  the  polyhedron  can  move  in  a  space-filling  way. 
A  consequence  of  this  result  is  that  the  structure  of  the  configuration  space  is  quite 
complex,  which  makes  finding  finger  tracking  algorithms  non-trivial. 

•  Finger  tracking  for  a  polyhedral  object.  Under  the  same  assumptions  as  above,  we 
have  shown  that  the  differential  motion  of  rhe  tracking  finger  is  given  by  a  4  x  4 
linear  system.  This  surprising  result  is  very  feasible  computationally,  especially  in  the 
context  of  simulation. 

•  Polygons  in  the  plane.  Our  newly  developed  framework  for  dexterous  manipulation 
has  been  used  to  express  earlier  results  [Rus90]  for  polygons;  these  results  were  origi- 
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naUy  obtained  in  a  more  ad  hoc  maimer.  Virtually  the  same  algorithm  used  to  deter¬ 
mine  the  configuration  space  for  polyhedra  can  be  used  to  determine  the  configuration 
space  for  polygons. 

•  Robustness  for  polygons.  The  planar  case  has  a  number  of  geometric  properties  that 
we  have  been  able  to  exploit  in  order  to  generate  rebus*  rotation  algorithms.  The 
uncertainty  inherent  in  the  real  world  makes  robustness  an  important  feature  for  any 
realistic  manipulation  algorithm.  Our  rotation  algorithms  are  robust  in  the  sense 
that  some  of  the  a  priori  knowledge  requirements  of  the  geometry  of  the  object  to 
be  rotated  and  the  necessity  for  precise  calculations  based  on  this  geometry  can  be 
replaced  by  sensing.  The  result  of  our  efforts  is  a  condition  on  the  geometry  of  the 
polygon  to  be  rotated  that  guarantees  robust  rotations  by  an  arbitrary  rotation  angle. 
We  ha,ve  shown  that  for  convex  polygons,  the  condition  can  be  checked  in  O(n)  time, 
with  O(nlogn)  preprocessing. 

We  are  currently  investigating  the  possibility  of  extending  the  robustness  results  from 
-the  planar,  polygon  case  to  the  3-dimensional,  polyhedron  case.  We  are  also  developing 
configuration  space  algorithms  for  3-dimensional  objects  with  curved  faces.  Another  area 
in  which  we  have  made  progress  is  the  experimental  verification  of  our  results.  We  are  using 
Newton,  the  simulator  for  rigid-body  dynamics  developed  by  our  group,  to  verify  our  finger 
tracking  algorithm  for  rotating  polyhedra. 
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A  Framework  for  the  Control  of 
Complex  Mechanical  Systems 


Dinesh  K.  Pai* 

Department  of  Computer  Science 
Cornell  University 
Ithaca,  NY  14853 


Abstract 

We  describe  an  approach  to  control  in  which  control  ac¬ 
tions  are  specified  as  weakly  as  desired.  We  use  large 
time-varying  sets  of  non-zero  measure  as  desired  goals  in¬ 
stead  of  specific  trajectories,  maintaining  that  we  do  not 
care  where  in  such  a  region  we  actually  are.  Inequality 
constraints  and  their  conjunctions  are  used  to  describe 
such  regions.  The  constraints  are  satisfied  at  run  time  to 
produce  the  control.  The  approach  has  been  successfully 
used  to  produce  human-like  walking  in  simulation. 

We  also  describe  an  implemented  programming  envi¬ 
ronment  for  this  approach.  We  discuss  the  representation 
of  control  computations  using  computational  graphs  and 
automatic  differentiation  for  efficiency.  Constraint  satis¬ 
faction  is  performed  using  a  fast  relaxation  method. 
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Figure  1:  Biped  walking  machine 


1  Introduction 

We  are  concerned  with  controlling  high  degree  of  free¬ 
dom  mechanical  systems  which  have  to  accomplish  several 
simultaneous  tasks.  Such  systems  include  robot  arms, 
multi-fingered  hands,  walking  machines,  mobile  robots, 
and  simulated  mechanical  systems  used  in  computer  an¬ 
imation.  The  system  typically  has  redundant  degrees  of 
freedom  for  each  task,  but  may  have  to  accomplish  a  large 
number  of  tasks  simultaneously.  The  system  may  also 
be  autonomous  and  reactive,  which  means  that  a  large 
amount  of  the  “programming”  will  be  done  at  execution 
time. 

The  complexity  of  the  system  implies  the  following: 

1.  The  ease  with  which  complex  tasks  are  expressed 
and  composed  is  critical. 

2  The  efficiency  of  the  control  computations  is  very 
important. 

3.  Simulation  of  the  mechanical  system  is  necessary  to 
gain  insight  into  the  control  programs  and  to  aid 
their  development. 

These  demands  are  frequently  contradictory.  We  propose 
a  framework  called  “Least  Constraint”  (abbreviated  as 

'Supported in  part  by  ONR  Grant  N 00014- &8K-0591,  ONR 
Grant  N00014-89J-1946,and  NSF  Grant  DMC- 86-17355.  The 
author  would  like  to  thank  L.-M.  Reissell  for  innumerable  dis¬ 
cussions  and  J  Cremer  for  comments  at  a  short  notice. 


LC)  which  we  believe  facilitates  the  expression  of  con¬ 
trol  programs  for  complex  mechanical  systems  and  is  ef¬ 
ficiently  implementable  as  well.  An  early  version  of  the 
framework  was  presented  in  [29]. 

Section  2  describes  the  LC  approach  and  how  motions 
are  expressed  in  it.  Section  3  discusses  the  data  struc¬ 
tures  for  representing  control  computations  in  LC  and 
the  algorithms  used  to  perform  these  computations  effi¬ 
ciently.  Section  4  describes  the  solution  of  the  inequality 
constraints  that  arise. 

2  The  Least  Constraint  Ap¬ 
proach 

2.1  Motivation 

Consider  the  problem  of  controlling  the  human-like  walk¬ 
ing  machine  of  Figure  1  to  walk  dynamically  in  three 
dimensions. 

One  approach  to  programming  such  a  task  is  to  pick 
some  periodic  trajectory  for  the  joints,  and  attempt  to 
track  it.  However,  it  is  not  clear  that  this  is  the  natural 
characterization  of  the  task.  Indeed  in  problems  of  this 
type,  a  major  gorJ  of  the  process  of  developing  a  control 
program  is  to  discover  the  task  requirements. 

We  would  instead  like  to  program  such  machines  in¬ 
crementally,  by  specifying  assertions  about  its  behavior. 
We  can  specify  several  requirements  for  walking,  for  in- 
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stance,  (i)  the  foot  should  clear  the  ground  during  the 
swing  phase  of  the  leg,  (ii)  the  swing  foot  should  be  moved 
to  a  location  suitable  for  dynamic  balance  by  foot  place¬ 
ment,  and  so  on. 

When  the  machine  is  controlled  to  satisfy  these  require¬ 
ments.  it  may  turn  out  that  the  requirements  were  inade¬ 
quate  —  for  example,  one  may  find  that  there  is  nothing 
to  prevent  knee  flexion  from  becoming  so  large  that  walk¬ 
ing  is  impossible.  In  this  case  one  would  like  to  modify 
the  existing  program  by  merely  adding  new  assertions:  for 
example,  by  adding  the  assertion  that  the  pelvis  should 
be  above  a  certain  height.  This  is  not  possible  in  current 
robot  programming  languages.  The  LC  framework  was 
designed  to  address  these  problems. 


2.2  Least  Constraint  Framework 


In  the  LC  framework,  motions  are  expressed  by  means  of 
time-  and  state-dependent  assertions.  These  assertions 
are  defined  using  inequality  constraints  which  describe 
the  set  of  allowed  states  as  a  function  of  lime.  The  con¬ 
straints  are  solved  at  run  time  to  produce  a  motion  sat¬ 
isfying  them. 

Since  complex  mechanical  systems  have  large  state 
spaces,  it  is  not  convenient  or  natural  to  express  all  of 
the  constraints  in  a  single  space.  For  instance,  the  walk¬ 
ing  machine  in  Figure  1  has  a  28-dimensional  state  space. 
For  convenience  of  expression,  users  define  derived  vari¬ 
ables  in  terms  of  the  basic  (e.g..  state)  variables  —  an 
example  of  this  is  the  definition  of  task  and  end-effector 
coordinates  for  robot  manipulators.  LC  generalizes  such 
constructions  to  allow  the  creation  of  arbitrary,  user  de¬ 
finable  quantities  which  are  natural  to  the  tasks  and  the 
constraints  being  expressed-  One  can  isolate  smail  groups 
of  variables  into  domains  on  which  to  focus.  For  example, 
the  foot  collision  constraints  in  the  above  walking  exam¬ 
ple  are  best  expressed  in  a  separate  foot  position  domain. 

In  LC.  users  define  a  domain  system,  {V,  .  i  6  /}, 
related  by  linking  functions 

Uj  :  Vi  —  Vj,  («,  j)  €  L  C  /  x  /,  (1) 

which  satisfy  the  basic  consistency  condition  that  all  di¬ 
agrams  of  the  following  form  commute. 

Iht  Imj 


Vi 


All  domains  Vi  are  connected  to  a  basic  domain  Vo  by 
compositions  of  linking  functions  : 


V0  if  ...  -  Vi . 


The  topological  properties  of  the  domains  V,  are  left 
somewhat  open  —  in  theory,  a  domain  system  could  be 
expressed  in  terms  of  arbitrary  manifolds;  in  practice, 
however,  the  domains  wiil  be  copies  of  7?".  for  varying 
n.  Non-Euclidean  domains,  such  as  SOi3y.  are  currently 
treated  using  their  coordinate  charts. 

A  motion  specification  in  LC  now  consists  of  a  system 
of  time-varying  inequality  constraints  Fa,  a  €  A.  on  the 
domains  V,.  here  each  constraint  is  expressed  by 

Pa=fa(rU).t)<0.  -2. 

where 

fQ  :  V,  x  H  —  ft.  a  £  A 

is  a  smooth  map,  and  r(t)  denotes  a  time-dependent  tra¬ 
jectory  in  Vi. 

The  interpretation  of  the  constraint  function  is  that 
the  system  is  controlled  to  make  the  specified  expressions 
PQ  ;=  /0(x(f).  t)  <  0  true  at  all  limes  t. 

Such  Pa  and  'heir  conjunctions 

P.=  f\Pa  >3, 

O 

are  executable  LC  motion  programs. 

As  an  example,  Figure  1  depicts  some  constraints  on 
the  position  of  the  foot  for  the  walking  machine. 

The  solution  of  the  constraints  is  achieved  by  produc¬ 
ing  a  trajectory  x(r)  6  Vo  such  that  the  derived  con¬ 
straints 


Pa  :=/«((/*.  o  —  o/0j)(x(t)).0  <  0 

obtained  by  lifting  t  le  original  constraints  using  ihe  link¬ 
ing  maps  are  satisfied  at  all  times  t.  In  LC.  this  is  done  at 
run  time  at  discrete  time  steps  at  every  time  step  r«. 
a  feasible  point  r(t»)  is  produced,  and  is  used  to  compute 
the  control  u. 

We  have  implemented  a  programming  environment  us¬ 
ing  Common  Lisp  and  the  X  window  system  to  develop 
and  execute  LC  programs.  In  this  environment,  users  de¬ 
fine  variables,  domains  and  linking  functions,  and  express 
constraints  in  these  domains.  Exact  partial  derivatives 
can  be  efficiently  computed  using  automatic  differentia¬ 
tion.  The  environment  interacts  with  a  simulator  of  rigid 
body  mechanics,  “Newton"  [8.14].  LC  accepts  models 
of  multi-body  mechanical  systems  described  in  Newton 
and  generates  the  computational  graph  (see  Section  3) 
for  the  simultaneous  computation  of  forward  kinematics 
of  user  specified  features  on  each  body,  differential  kine¬ 
matics,  and  inverse  dynamics.  The  control  programs  c\n 
be  tested  by  sending  the  control  output  produced  by  LC 
to  the  simulator.  Tools  are  provided  for  debugging  con¬ 
trol  programs  and  for  visualizing  constraints.  Finally, 
the  control  programs  can  be  translated  into  straight  line 
programs  in  other  languages  (currently  Lisp  and  C)  for 
efficiency. 


Briefly,  the  motivation  for  using  domain  systems  is  that 
they  allow  a  constraint  on  a  subdomain  Vi  to  be  lifted 
to  an  equivalent  constraint  on  the  basic  domain  Po  using 
compositions  of  Unking  functions. 


2.3  Examples 

We  have  deUberately  left  open  issues  such  as  how  the  sys¬ 
tem  is  used,  the  nature  of  the  constraints,  and  the  choice 


Figure  2:  Toleranced  move  near  singularity. 


of  basic  variables.  LC  is  intended  to  be  a  general  pro¬ 
gramming  framework,  rather  than  a  solution  to  a  specific 
problem.  In  this  section,  we  present  examples  of  the  ap¬ 
plication  of  LC  to  robot  programming.  The  examples  are 
simple  for  the  sake  of  presentation,  but  illustrate  the  use 
of  LC. 

Example  1.  Figure  2  depicts  a  2-link  planar  robot  ma¬ 
nipulator  which  is  to  be  controlled  to  move  along  a 
straight  line  passing  through  the  singular  configuration 
(corresponding  to  the  end-effector  location  (0.0)),  from 
the  point  (1.0, 0.0)  to  (-1. 0,0.0)  with  an  average  speed 
of  at  least  v.  We  assume  that  this  is  a  pure  kinematic 
problem  —  i.e..  the  basic  variables  are  the  joint  angles 
1 1  and  tj,  which  are  tracked  by  a  separate  low  level  con¬ 
troller.  We  are  allowed  a  tolerance  of  (  around  the  nom¬ 
inal  straight  line. 

The  assertion  that  the  end-effector  position  (rj.yj)  is 
within  <  of  the  nominal  straight  line  is  expressed  by  the 
constraint  functions 

f'op(li)  =  yi-e  (4) 

fbot(Vi)  =  -Vi-i  (5) 

and  finally  the  end-effector  is  ‘‘pushed”  along  the  tube  by 
specifying  the  constraint  function 

/pu>/i«r(X2,0  =  (*2  -  1.0)  -  tit  (6) 

Our  constraint  satisfaction  procedure  is  robust  near  sin¬ 
gularities  (Section  4)  anci  the  robot  executes  the  pose¬ 
changing  motion  shown. 

Example  2.  Figure  3  depicts  the  task  of  placing  an  ob¬ 
ject  in  the  robot’s  hand  on  the  table.  In  a  robot  pro¬ 
gramming  language  such  as  VAL  (34],  this  will  have  to 
be  expressed  by  arbitrarily  commanding  a  particular  mo¬ 
tion  to  an  arbitrary  point  on  the  table.  In  LC,  the  task 
is  programmed  r,j  placing  constraints  in  the  hand  posi¬ 
tion  domain.  The  hand  is  constrained  to  stay  within  the 
boundaries  of  the  table  by  the  constraints  fi  and  ft.  The 
constraint  ft  moves  down  at  constant  speed  v  from  height 
h.  so  that  af*e;  a  time  of  h/v  the  end-effector  is  on  the 
table  top.  Note  that  the  exact  location  on  the  table  top 
is  unspecified. 


If  we  now  wish  to  place  the  object  at  a  particular  lo¬ 
cation  on  the  table,  this  is  achieved  by  adding  the  cone- 
shaped  constraint  marked  f\  in  Figure  3.  If  it  is  required 
that  the  center  of  mass  of  the  arm  be  above  a  small  re¬ 
gion  of  the  base  during  the  motion,  a  new  variable  r,-, 
is  defined  in  terms  of  the  joint  angles  and  base  position, 
and  the  constraints  ti  <  <  ir. 

Example  3.  We  have  used  the  LC  framework  to  success¬ 
fully  program  the  human-like  machine  in  Figure  1  to  wajk 
dynamically  in  three  dimensions  [30j.  In  each  qualitative 
state  of  the  machine,  such  as  "standing  on  the  left  leg 
while  stepping  with  the  right.”  constraints  are  imposed 
to  achieve  a  large  number  of  simultaneous  tasks  such  as 
foot  placement  for  dynamic  balance,  torso  orientation, 
maintenance  of  pelvis  height,  collision  avoidance  with  the 
ground  during  swing,  inter-leg  collision  avoidance,  joint 
limit  avoidance,  etc. 

2.4  Discussion 

One  of  the  main  advantages  of  the  LC  framework  is  that  it 
enables  and  perhaps  encourages  subtasks  to  be  expressed 
weakly.  This  reduces  the  number  of  arbitrary  decisions 
which  have  to  be  made,  such  as  arbitrarily  picking  a  tra¬ 
jectory  in  Example  2.  The  program  reflects  the  user's 
intention  better  and  is  easier  to  maintain.  It  is  also  eas¬ 
ier  to  do  “redundancy  maintenance."  i.e..  to  retain  the 
excess  degrees  of  freedom  available  to  perform  a  task. 

LC  introduces  certain  object-oriented  features  to  robot 
programming.  Motion  programs  are  easy  to  combine,  in¬ 
herit,  and  specialize  for  new  tasks.  This  is  particularly 
important  for  complex  systems  for  which  programs  are 
developed  incrementally  (see  Examples  2  and  3.) 

Typical  users  will  use  libraries  of  higher-level  com¬ 
pound  programs  for  common  LC  idioms.  The  following 
is  a  brief  list  of  such  idioms. 

e  The  obttacle-.  These  constraints  encode  the  free 
space  in  the  environment  te.g.,  [21.4.10,7]).  Care 
must  be  taken  to  grow  the  constraints  or  to  impose 
velocity  constraints  in  order  to  account  for  the  brak¬ 
ing  characteristics  of  the  manipulator  near  the  ob¬ 
stacle. 
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•  The  interval:  This  is  a  special  case  of  the  obstacle 
idiom,  and  restricts  the  range  of  a  variable,  as  in  the 
case  of  joint  angle  limits  [17],  speed  limits,  etc. 

•  The  pusher.  This  a  time-dependent  constraint  which 
moves  the  system  in  a  given  direction,  without  re¬ 
stricting  motion  orthogonal  to  this  direction.  Equa¬ 
tion  6  is  an  example. 

•  The  funnel  [23]:  Here  the  constraints  define  a  set 
which  contracts  over  time.  Thus  the  system  can  be 
brought  to  a  desired  configuration  without  overly 
restricting  its  trajectory.  A  canonical  example  of 
a  funnel  is  a  contracting  ball. 

•  The  toleranced  move :  A  moving  ball  with  fixed  ra¬ 
dius  is  an  example.  Example  1  provides  another  in¬ 
stance  of  a  toleranced  motion. 

2.5  Relationship  to  other  approaches 

Current  approaches  to  programming  complex  mechanical 
systems  may  be  broadly  classified  as  follows: 

•  Explicit  approaches:  These  consist  of  approaches 
which  allow  users  to  explicitly  specify  the  mo¬ 
tion  of  the  robot,  e.g..  by  prescribing  trajectories. 
These  include  robot-level  programming  languages 
(e  g  ,  [34,33]).  While  these  approaches  allow  fine 
control  of  the  motion,  they  are  hard  to  program  and 
force  users  to  make  arbitrary  decisions  in  order  to 
execute  a  motion. 

•  Implicit  approaches:  In  these  approaches,  the  users 
only  specify  the  high-level,  global  goals  of  the  mo¬ 
tion  and  the  system  plans  a  motion  that  achieves  the 
goals.  These  include  the  approaches  of  motion  plan¬ 
ning  (e.g.,  [21. 15.22. 11.S.20])  and  optimal  trajectory 
planning  [3.32.6].  These  methods  are  powerful  when 
they  are  well  matched  to  the  problem.  However  this 
is  frequently  not  the  case;  for  example,  it  may  not 
be  important  that  the  trajectory  minimize  a  spe¬ 
cific  functional  such  as  energy  along  the  trajectory. 
More  generally,  these  approaches  do  not  facilitate 
modification  of  the  planned  motions  by  the  users. 
Finally,  they  are  typically  computationally  expen¬ 
sive  and  need  to  be  executed  off-line. 

LC  is  an  intermediate  approach  between  the  explicit 
and  implicit,  sharing  some  features  of  both. 

LC  is  higher  level  than  the  explicit  approaches  — 
one  need  not  specify  motions  explicitly  but  rather  more 
weakly  as  a  set  of  constraints.  On  the  other  hand,  by 
moving  the  constraints  and  restricting  the  feasible  set. 
one  has  a  degree  of  explicit  user  control  on  the  motion. 

LC  is  lower  level  than  the  implicit  approaches.  It 
has  no  built-in  application  specific  knowledge.  The  con¬ 
straints  are  satisfied  locally,  and  LC  cannot  guarantee 
that  the  motion  generated  at  one  time  will  not  cause  a 
failure  (e  g  .  there  may  be  no  feasible  point  in  a  small 
neighborhood  of  the  current  state)  at  some  time  in  the 
future  An  additional  planning  iayer  may  be  necessary  to 
avoid  such  situations.  On  the  other  hand,  for  problems 
such  as  controlling  reactive,  autonomous  systems  (e.g.. 
[5]),  this  lack  of  guarantees  is  not  a  critical  issue. 


There  is  a  close  connection  between  satisfying  con¬ 
straints  and  avoiding  obstacles  --  obstacles  are  ’hard' 
inequality  constraints.  Conversely,  one  can  think  of  con¬ 
straints  as  being  virtual  obstacles  m  abstract  spaces, 
which  can  change  and  move  over  time  under  the  pro¬ 
grammer’s  control.  In  particular,  our  approach  shares 
several  features  with  the  use  of  artificial  potential  fields, 
proposed  oy  Khatib[!7]  (see  also,  for  example.  fl2.2T.lV 
However,  there  are  important  differences. 

First.  LC  generalizes  the  notion  of  obstacles  to  con¬ 
straints  in  arbitrary,  user-defined  domains.  This  is  sup¬ 
ported  by  a  programming  framework  to  describe  these 
constraints  conveniently,  and  a  constraint  satisfaction 
system  which  solves  the  constraints.  Thus  LC  should  be 
applicable  to  a  broad  range  of  motion  control  tasks. 

Another  difference  is  the  specification  of  motions  out¬ 
side  the  natural  constraints  imposed  by  obstacles  and 
joint  limits.  In  potential  function  approaches  the  mo¬ 
tion  is  specified  by  constructing  a  scalar  function  o  such 
that  the  system  behaves  like  a  gradient  dynamical  system 
with  o  as  potential.  The  specification  of  a  motion  using 
a  potential  function  is  concise  —  a  single  scalar  function 
encodes  global  dynamic  behavior  of  the  system.  How¬ 
ever.  this  has  the  drawback  that  it  is  difficult  for  users 
to  specify  potential  functions  for  complicated  behaviors. 
Thus  the  potential  functions  encountered  in  the  literature 
are  extremely  simple  or  are  generated  by  special  purpose 
planning  programs  (see  [19]).  In  LC  motions  are  specified 
by  time  varying  constraints.  Each  constraint  has  an  intu¬ 
itive  meaning  as  an  assertion.  The  construction  of  com¬ 
plicated  potential  functions  by  adding  simpler  potentials 
together  does  not  necessardy  result  in  easily  predictable 
behavior:  but,  joining  two  constraints  will  always  pro¬ 
duce  motion  which  satisfies  both  constraints.  The  state 
is  manipulated  by  time-varying  constraints  in  a  manner 
reminiscent  of  pushing  operations  [24]  and  may  be  viewed 
as  a  generalization  of  pushing  to  user-defined  spaces. 

We  believe  that  these  features  make  constraints  easier 
to  specify  and  visualize  than  potential  functions. 

?  Representing  Computations 

3.1  Computational  Graph 

The  domain  system  of  Section  2.2  is  implemented  in  LC  as 
a  computational  graph  (or  Kantorovich  graph)  [31.16}-  In 
a  typical  constraint  satisfaction  computation,  one  needs 
to  compute  the  value  of  each  constraint,  and  the  gradi¬ 
ents  of  the  violated  constraint  functions.  These  compu¬ 
tations  are  efficiently  performed  using  the  technique  of 
automatic  differentiation.  Computations  of  derived  vari¬ 
ables,  c.pecially  for  kinematics,  differential  kinematics, 
and  dynamics,  contain  many  common  subexpressions,  the 
elimination  of  which  is  also  a  major  source  of  efficiency. 
The  computational  graph  is  a  useful  data  structure  for 
achieving  these  goals. 

A  computational  graph  can  be  described  as  follows  isec 
Figure  4  for  example).  A  function  /  .  —  ft*  is  said  to 

be  factorable  if  every  component  of  /  is  a  function  com¬ 
putable  from  the  basic  and  derived  variables  by  means  of 


i/2 


Figure  4:  Computational  graph  for  2-link  robot 


a  finite  sequence  of  elementary  operations. 

We  represent  a  factorable  function  as  a  graph  as  fol¬ 
lows:  let  V[  be  the  set  of  m  vertices  corresponding  to 
the  m  coordinates  of  the  domain  of  /.  k'o  be  vertices 
corresponding  to  the  n  coordinates  of  the  range,  and 
Vt-  be  the  intermediate  quantities  in  the  computation 
of  /.  Let  G(V.  E)  be  a  directed  graph  with  vertex  set 
V'  =  {  V/  uFyU  Vo)  and  edge  set  E.  Let  i>  be  the  set  of 
basic  operations.  With  each  vertex  we  associate  an  oper¬ 
ation  by  the  function  .  V  —  v  with  the  understanding 
that  the  operation  ui(c)  is  applied  to  the  ordered  list  of 
immediate  predecessors  of  u  in  G. 

Figure  4  depicts  a  computational  graph  for  the  y- 
coordinate  of  the  end-effector  position  of  a  2-link  robot 
manipulator. 

Vertices  of  the  graph  are  grouped  together  to  form  do¬ 
mains,  and  the  subgraph  connecting  two  domains  defines 
their  linking  function.  The  basic  domain  Vo  consists  of 
all  the  independent  variables  in  the  system's  computa¬ 
tional  graph,  i.e.  the  set  of  all  vertices  with  in-degree  0, 
which  are  not  marked  as  constants.  Frequently,  the  basic 
domain  is  identical  to  the  state  space  of  the  system. 


i/2  dt7 


Figure  5:  Augmented  computational  graph  for  2-lins 
robot 

ample  in  Figure  4.  In  the  example  dyo/dt  i  is  the  same 
as  zj,  the  r-coordinate  of  the  end-effector  position,  and 
need  not  be  recomputed. 

In  the  reverse  mode,  the  additional  cost  of  computing 
all  the  partials  of  a  function  is  very  small.  An  upper 
bound  was  derived  by  Baur  and  Strassen  (2.251. 

Let  /  be  a  rational  function  of  m  variables  z; . zn 

for  which  the  computational  graph  has  iqG)  vertices,  of 
which  p(G)  are  muitiplications/divisions.  Then 

r(G')  <  5r(G) 

11(C)  <  3  p(G)  (7) 

Thus  the  cost  of  computing  all  the  partials  of  a  scalar 
function  is  at  most  a  small  constant  multiple  of  the  cost 
of  merely  evaluating  the  function.  The  extension  to  tran¬ 
scendental  functions  and  exponentiation  is  straightfor¬ 
ward  [25].  Our  experience  is  that  the  above  bounds  are 
conservative.  Also,  it  is  frequently  the  case  that  the  func¬ 
tion  has  been  computed  before  the  gradient  is  required,  in 
which  case  the  cost  of  computing  the  gradient  is  reduced 
by  r(G)  (or  p(G)). 


3.2  Automatic  Differentiation 

Automatic  differentiation  is  a  technique  for  efficiently 
computing  exact  derivatives  of  factorable  functions  (see 
[31,16,13]  for  recent  surveys).  It  differs  from  standard 
symbolic  differentiation  in  that  by  using  the  underlying 
computational  graph  data  structure,  the  growth  of  com¬ 
mon  subexpressions  due  to  differentiation  is  automati¬ 
cally  controlled;  it  differs  from  numerical  differentiation 
in  that  the  method  is  exact  and  results  in  no  loss  of  sig¬ 
nificant  digits. 

We  have  implemented  both  forward  and  reverse  modes 
of  automatic  differentiation.  The  computation  of  par¬ 
tial  derivatives  is  implemented  as  augmentation  of  the 
original  computational  graph  to  produce  a  new  graph 
G'  that  computes  both  the  function  and  its  derivatives. 
This  leads  to  additional  savings  since  the  expressions  for 
the  partials  are  frequently  already  present  in  the  com¬ 
putational  graph  and  are  found  using  simple  expression 
matching  algorithms.  Figure  5  depicts  the  augmented 
computational  graph  of  the  2-link  robot  kinematics  ex¬ 


4  Constraint  Satisfaction 

LC  separates  the  specification  of  constraints  and  the  tech¬ 
nique  used  for  satisfying  them.  Thus  different  constraint 
satisfaction  algorithms  can  be  used.  We  describe  here  a 
non-linear  extension  of  the  relaxation  method  [it.  which 
can  efficiently  handle  constraints  in  multiple  domains 
The  method  is  local  and  iterative.  The  local  nature  of 
the  method  makes  it  fast  enough  to  be  implementable  in 
real  time.  It  exploits  the  knowledge  of  a  good  starting 
point  for  the  iteration  —  the  solution  to  the  constraint 
satisfaction  problem  at  the  previous  time  step. 

4.1  Constraints  in  a  single  domain 

First  co.isider  the  case  in  which  all  constraints  ate  ex¬ 
pressed  in  one  domain.  The  method  we  use  for  consttamt 
satisfaction  is  similar  to  the  relaxation  method  for  linear 
inequalities  [1],  but  is  extended  to  nonlinear  inequalities 
using  Newton’s  method. 
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Let  V  =  Sm  be  the  domain,  z  S  2>,  and  let  the  system 
of  inequality  constraints  be 

/,(r)  <0,i=  ! . n.  (8) 

Let  zk  oe  the  Jfcth  iterate  at  time-step  l,  and 

-  ■  M 

-  - 

Then  each  iteration  is: 

i  For  each  i  =  I _ n  compute  /,(z*).  If  aii  /,(z*)  < 

0.  terminate. 

2.  Else,  let  j  be  such  that  ds  =  max,  d;. 

**'rl  =  r*  -  P[d3  +o)n.  (11) 

p  6  (0, 2]  is  the  relaxation  parameter,  and  o  is  a 
small  offset. 

The  behavior  of  the  above  algorithm  has  been  analyzed 
in  the  case  of  linear  /;  by  [1.26}.  The  algorithm  s  then 
globally  convergent,  with  positive  offsets  o  having  the  ef¬ 
fect  of  inducing  convergence  in  a  finite  number  of  steps. 
The  behavior  in  the  non-linear  case  is  similar  in  a  suffi¬ 
ciently  small  neighborhood  of  the  feasible  set.  The  global 
convergence  property  is  lost,  but  this  is  not  a  serious  prob¬ 
lem  since  we  are  typically  solving  the  inequalities  starting 
near  the  feasible  set.  Finally,  the  convergence  can  be  lin¬ 
ear  with  a  large  convergence  ratio  when  the  inequalities 
are  ill  conditioned-  Chocring  relaxation  parameter  p  >  1 
can  significantly  improve  convergence.  We  are  currently 
investigating  the  tradeoffs  with  other  methods  that  have 
super-linear  convergence  but  are  more  expensive  per  it¬ 
eration. 

The  algorithm  works  well  in  practice.  Step  1  of  each 
iteration  can  be  computed  in  parallel  to  check  violated 
constraints.  The  gradient  required  in  equations  9  and 
10  is  efficiently  computed  using  automatic  differentiation, 
and  utilizes  the  effort  already  expended  in  computing  /;. 

4.2  Constraints  in  multiple  domains 

Since  LC  users  specify  constraints  in  different  domains, 
possibly  related  by  linking  functions,  the  single  domain 
case  above  has  to  be  extended  to  handle  multiple  do¬ 
mains  Consider  the  case  of  two  domains  V\  and  Vz.  and 
further  assume  that  V\  is  defined  in  terms  of  Vz  by  a 
smooth  linking  map  l :  Th  —  V\ .  Let  fl  be  a  constraint 
in  V i  and  /s  be  a  constraint  in  Pj. 

It  is  computationally  advantageous  to  perform  the  de¬ 
scent  in  Vz-  We  lift  /l  to  as  equivalent  constraint 
P  =  fl  of  in  Vz-  Thus  V/5  =  D/rV/1.  making  the 
method  robust  sear  singularities. 

Also,  it  is  not  necessary  to  explicitly  compute  the  ma¬ 
trix  DlT  —  we  directly  compute  V/1  by  differentiat¬ 
ing  P  using  automatic  differentiation  saving  the  cost  of 
the  dim  V\  xHim  Vz  multiplications  for  the  matrix- vector 
product  DlT^p. 


This  method  is  preferable  to  methods  which  choose 
the  steepest  descent  and  the  estimate  of  the  distance  :n 
the  specification  domain  V\.  Although  the  latter  choice 
makes  it  easier  for  the  user  to  anticipate  the  behavior 
of  the  constraint  satisfaction  algorithm,  it  is  also  mote 
expensive,  in  effect  requiring  the  computation  of  l~‘  or 
Dl~‘ .  This  also  leads  to  numerical  problems  when  the 
Jacobian  matrix  Dl  is  singular  or  til  conditioned.  This  is 
important  in  robotics  applications,  where  motion  in  the 
vicinity  of  kinematic  singularities  is  sometimes  desirable 
or  inevitable  [28]. 

The  general  situation  involves  constraints  expressed  in 
several  related,  and  possibly  over.appmg.  domains  and  is 
handled  similarly  by  lifting  each  constraint  function  /,  ;n 
domain  V,  to  the  basic  domain  Vz  by  using  a  composition 
of  the  linking  functions  provided. 

4.3  Cost  per  iteration 

Let  c,  be  the  cost,  in  mu.tiplications.  of  computing  each 
of  the  n  constraint  functions  /,.  The  cost  of  computing 
all  the  /,  is  V  <  c,.  where  the  *<"  case  arises  if  the 

/■  share  computations.  The  test  of  constraint  violation 
in  Step  I  of  the  iteration  costs  V  and  this  is  fixed  by 
specification  of  the  constraints  by  the  user.  Let  A  = 
{:’!/,  >  0}  be  the  active  or  violated  constraints.  Then  the 
update  in  Step  2  costs 

L<^2 c.  (121 

•eA 

Here  we  assume  teat  c,  >  dim  Vo  and  terms  of  the  order 
of  dim  VQ  are  ignored.  The  cost  of  computing  the  gradient 
is  actually  bounded  by  3c,.  but  of  this  c,  has  already  bees 
accounted  for  :n  V*.  The  total  cost  of  each  iteration  -s 
V~U. 

Note  that  U  is  very  small  when  only  a  few  constraints 
are  violated,  and  is  at  most  twice  the  cost  of  merely  check¬ 
ing  if  any  constraints  are  violated. 

5  Conclusions 

The  basic  idea  of  the  approach  to  control  described  m 
this  paper  is  'Least  Constraint^  —  by  which  we  mean 
that  we  specify  control  actions  as  weakly  as  desired.  This 
permits  motions  specifications  to  better  reflect  the  user  s 
intentions,  and  makes  the  programs  easier  to  maintain. 
Another  contribution  of  this  work  is  our  use  of  computa¬ 
tional  graphs  for  the  efficient  computation  of  a  constraint 
function  and  it's  derivatives.  This  leads  to  an  inexpensive 
constraint  satisfaction  algorithm  that  has  the  auxiliary 
benefit  that  it  is  robust  near  singularities. 

We  have  used  the  LC  framework  on  a  small  range 
of  problems,  including  the  task  of  programming  a  sim¬ 
ulated  human-like  biped  machine  to  walk  dynamicaly 
Out  experience  indicates  that  LC  is  a  useful  programming 
approach  for  complex  systems  with  several  constraints 
However,  the  method  can  perform  pooriy  when  the  con¬ 
straints  are  ill-conditioned  and  the  feasible  set  a  small 
and  disconnected.  For  such  cases,  we  are  considering 


methods  with  super-linear  convergence  properties  and  for 
selecting  viable  connected  components  of  the  feasible  set 
using  simulation. 

Finally,  we  would  like  to  suggest  that  LC  may  be  a 
suitable  language  in  which  high  level  planners  express  mo¬ 
tions.  The  plans  produced  can  then  be  easily  augmented 
by  sensors  or  tweaked  by  users  if  nec.-ssaiy. 
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Abstract 

Recently,  implicit  patches  have  emerged  as  alternative 
modeling  primitives  for  three  dimensional  objects  In 
designing  three  dimensional  models,  one  often  encoun¬ 
ters  various  shape  requests.  This  paper  develops  tech¬ 
niques  for  satisfying  such  requests  through  shape  con¬ 
trol.  Jn  particular,  we  show  how  to  achieve  the  convexity 
of  quadric  patches  or  cubic  patches. 

1  Introduction 

The  end  goal  of  geometric  modeling  is  to  design  and 
to  manipulate  three  dimensional  models  represented  by 
free-form  surfaces.  Traditionally,  free-form  surfaces  are 
built  from  parametric  patches.  Parametric  patches  are 
successful  as  far  as  design  and  rendering  are  considered, 
but  manipulating  three  dimensional  models  with  para¬ 
metric  pathes  poses  fundamental  difficulties.  For  exam¬ 
ple,  parametric  patches  are  not  closed  under  sweeping 
and  convolution.  The  intersection  of  two  parametric 
patches  are  extremely  difficult  to  represent  and  evalu¬ 
ate  [HK86]. 

One  way  to  avoid  these  problems  is  to  build  fiee- 
form  surfaces  from  low-degree  implicit  patches.  Implicit 
patches  are  closed  under  all  common  operations  in  ge¬ 
ometric  modeling  [Baj88],  and  the  intersections  of  low- 
degree  implicit  patches  can  be  computed  efficiently[OSS]. 
Recent  research  shows  that  quadric  and  cubic  implicit 
patches  are  flexible  enough  for  building  arbitrary  tluee 
dimensional  models  [Guo90,  Guo91). 

A  major  reason  that  parametric  patches  have  be¬ 
come  so  popular  in  computer  graphics  is  their  good 
shape  control  properties.  In  this  paper,  we  tachle  the 
shape  control  issues  of  implicit  patches.  Using  Bernstein- 
Bezier  representation  of  polynomials,  we  can  control  the 
shapes  of  implicit  patches  through  manipulating  theii 
control  points. 

"This  work  is  supported  by  DARPA  under  ONR  contract 
N00014-86K-0591,  NSF  Grant  DMC-8G- 17355.  and  ONR  Grant 
N00014-89J-  1946. 


In  designing  free-form  surfaces,  one  often  encounters 
various  shape  requirements,  such  as  a  nice  pattern  of  re¬ 
flection  lines  and  lestrictious  on  the  minimum  radius  of 
Cnivatuic.  Among  nil  the  shape  leqniiements.  convex¬ 
ity  is  thv  most  basic  and  the  most  fieqnently  leq nested 
one.  In  this  papei.  we  show  how  to  manipulate  the  con- 
tiol  points  of  a  quadlie  patch  or  a  cubic  patch  so  that 
the  patch  become  convex. 

1.1  Previous  work 

I.ow-degtee  implicit  sui faces  ate  extensively  used  in  the 
existing  solid  modeling  and  giaphics  systems  as  mod¬ 
eling  piimitives  [RY83].  and  geometric  opetations  on 
low-degiee  implicit  surfaces  aie  well  understood  [OSS]. 
Implicit  surfaces  aie  also  very  useful  in  sitiface  fitting 
[PK89]  and  blending  [HOST.  MS85,  IIII87,  BliS'J], 

Many  authors  have  addressed  the  shape  contiol  of 
implicit  patches  [Sed«5,  VVMVV8(5,  BVV90).  In  part ic- 
ulai.  Bloomcuthal  and  Wyvill  [BVV'IO]  discussed  shape 
contiol  using  skeletons,  and  Sedeiberg  pointed  out  that 
the  Bei nstein-Beziei  representation  ate  suitable  forcon- 
t tolling  the  shapes  of  implicit  patches  [SedS.>] 

1.2  Overview 

This  papei  is  oigaui/td  as  follows.  Aftci  giving  some 
buchgioiiud  iufuim.it ion  in  Section  2,  we  dcsciibe  the 
basic  shape  contiol  tiihlliqlies  ill  Section  3.  Section  1 
shows  how  to  achh  u  tin  convexity  of  quadtic  patches 
and  cubic  patches 

2  Bernstein-Bezier  representation 

(liven  ;t  lelrahedion  \  with  veitices  X\ .  X?.  X3.  and  X4. 
one  can  t \|>ies>  any  point  p  in  >pace  <\* 
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Figure  1:  Cubic  control  points 


The  tuple  (ri,  12,73,74)  is  called  the  barycenttic  coor¬ 
dinate  of  p.  The  barycentric  coordinates  are  linearly 
related  to  Cartesian  coordinates,  so  any  implicit  poly¬ 
nomial  surface  may  be  expressed  in  barycentric  coot  di- 
nates  via  a  linear  change  of  variables. 

For  a  non-negative  integer  tuple  A  =  (Ai.Aj.As.A.)) 
with  |A|  =  Vi.  A;  =  n,  the  Bernstein  polynomial  for 
A  is 


Using  Bernstein  polynomials,  one  can  uniquely  repre¬ 
sent  any  polynomial  /  of  degree  <  n  as  follows. 


f(r)  =  Y, 

m =» 


The  b\’s  are  referred  to  as  the  control  points  of  the  poly¬ 
nomial  /  and  its  surface  S(f)  =  (x|/(x)  =  0}.  The  con¬ 
trol  points  of  a  cubic  polynomial  are  shown  in  Figure  1. 
The  following  lemma  is  very  useful. 

Lemma  1  If 


f(x )  =  Y 

W=* 


and  ckt,  =  0,  then 


3  The  basic  techniques  for  shape  con¬ 
trol 

An  implicit  patch  is  defined  as  the  zeto  conlout  of  a 
polynomial  f  inside  a  letrahcdion  [x|X2X.iX1]1. 

f(r)=  Y  bxBUr). 


where  r  is  the  baivceniric  coordinate  defined  by  the 
tetrahedron  [xiX2XjX.|].  The  basic  idea  of  shape  con- 
tiol  is  to  expiess  the  geometiic  propeities  the  implicit 
patch  in  teuns  of  the  control  points  so  that  one  can 
achieve  shape  objectives  by  lequiiing  the  contiol  points 
to  satisfy  certain  constraints. 

To  get  a  feel  foi  the  effects  of  the  contiol  points,  we 
study  a  univuiiatc  cubic  poly  \l  /. 

f(  It )  =  l>3 0  Bio  +  I>21  B'ii  I  2  Bio  +  boo  Be, 3 

The  value  of  the  function  /  over  [0 . 1  ]  and  the  convex  hull 
of  the  points  (0. 6ao  )*  (1/3. 621).  (2/3.  612),  and  (l.ios) 
aie  shown  in  Figme  2.  The  functions  #21-  $12- 
and  /Jq3  aie  shown  m  Figme  3. 

Fiom  these  liguies.  we  can  see  the  following. 

1.  At  the  end  points  of  the  interval  (O.lj.  the  control 
points  boo  and  boo  equal  to  the  function  values  of 
/■ 


c(k~l)e‘+e>  —  (^7/(x<)  ixJ—X‘)< 

for  j  =  1,2, 3, 4. 

Proof:  From  [Dah86]  , 

(x;  -Xj,V/(x)) 

=  k  Y^/  (CA+e‘  —  c\+e>  )B\  l(r) 
|A|=fc-l 

Letting  x  =  x, ,  we  prove  the  lemma.  £ 


2.  The  giadient  of  /  at  the  end  points  of  the  interval 
[0.1]  aie  determined  by  boo-  /'21  ■  b  12.  and  boo- 

3.  The  contiol  point  boo  has  a  effect  on  the  value  of 
f(u  1  foi  all  11  except  «  =  0.  and  the  effect  is  the 
st longest  near  11  =  I  Similai  statement  can  be 
made  about  olliei  contiol  points. 

All  these  lehrtions  between  the  contiol  points  and  the 
piopeilics  of  the  polynomial  /  geneiah/e  to  tnvariate 
polynomials 

1  \\V  <leiiotc  l>\  [X| . 


.  Xi  j  (lie  comex  hull  of  j  Xi .  .  Xj,  } 


•  control  point 

Figure  2:  Function  values  of  a  univaiiate  cubic  function 


Having  understood  the  elTect  of  the  control  points, 
we  use  the  control  points  to  control  the  shape  of  the 
surfaces.  Consider  the  problem  of  interpolating  points 
and  lines  in  space  by  a  surface  S(/).  Since  the  values 
of  /  at  a  vertex  of  the  tetrahedron  [X1X2X3X4]  is  equal 
to  the  value  the  control  point  at  the  vertex,  setting  the 
control  point  to  zero  forces  the  S(f)  to  pass  through 
the  vertex.  This  method  of  interpolating  a  points  can 
be  generalized  to  a  method  of  interpolating  the  edges  of 
the  tetrahedron  [X1X2X3X4]. 

Moving  on  to  the  problem  of  controlling  the  tangent 
plane  of  S(f),  we  consider  the  tangent  plane  of  S(f)  at 
the  vertex  Xi.  The  tangent  plane  at  Xi  is  defined  by  its 
the  gradient  V/(x  1).  From  Section  2,  we  know  that 

b(k-l)ei+c>  =  £(V/(x,),Xj  —  x, ),  j  =  2.3.4. 

Since  the  vectors  xj  -  Xi  (j  =  2,  3,4)  are  three  linearly 
independent  vectors,  the  above  relation  implies  that  the 
control  points  next  to  Xi  completely  determine  the  gra¬ 
dient  V/(x  i). 

More  sophisticated  examples  of  shape  control  arc- 
easy  to  come  by.  The  restriction  of  /  to  an  edge  of 
the  tetrahedron  (X1X2X3X4]  is  a  univariate  polynomial. 
If  the  control  points  on  the  edge  are  all  positive  or  all 
negative,  then  the  surface  S(f)  does  not  intersect  the 
edge.  Otherwise,  the  surface  S(f)  intersects  the  edge- 
exact  once  if  there  is  exactly  one  sign  change  111  the  list 
of  control  points  along  the  edge.  Similar  statemeius  ca  1 
be  made  for  the  fa"s  of  [xlX2X3X<]. 

4  Achieving  convexity 

As  an  application  of  the  techniques  for  shape  control, 
we  derive  the  convexity  condition  of  an  implicit  patch  in 
terms  of  its  control  points.  Throughout  this  section,  we 
concentrate  on  the  implicit  parch  defined  as  the  portion 
of  a  surface  S(f)  inside  a  tetrahedron  V  =  [X1X2X3X4]. 

The  reader  is  familiar  with  convex  objects  as  a  set 
of  points  in  three  dimensional  space  such  that  the  line 
segment  connecting  two  points  in  the  set  is  contained 
in  the  set.  Convex  surfaces  are  often  defined  as  surfac  -s 
whose  Gaussian  curvatures  are  positive  over  the  entile 


surface.  Convex  objects  and  convex  surfaces  are  related 
in  t Ii.it  if  a  convex  sulfate  is  closed  and  it  bounds  a 
point  set  vvith  finite  volume,  then  the  point  set  is  a 
convex  object. 

Defining  a  convex  sutface  in  terms  of  Gaussian  cur- 
vatuic  is  not  convenient  when  dealing  with  implicit  sur¬ 
faces  So  we  use  the  definition  of  convex  suifaces  in 
let  tils  of  the  tangent  planes.  Let  an  implicit  surface 
-'?(/)  have  a  tangent  plane  ■'i(Px)  at  point  x  G  S(f). 
The  mii face  .$(/)  is  convex  at  the  point  x  if  the  surface 
5(/)  *s  111  tlie  half  space  bounded  by  5(  P\ )  and  pointed 
to  bv  -V/(x).  An  implicit  patch  is  convex  if  its  pri- 
maiy  mii  face  is  convex  at  even  point  on  the  implicit 
patch. 

Notice  the  iclationship  between  the  convexity  of  the 
siii face  S(f)  and  the  convexity  of  the  polynomial  /,  A 
polynomial  /  is  convex  ovei  the  tetiahedion  1’  if  for  any 
two  points  x  and  y  in  the  (ciiahc-dton. 

)<  j(/(x)  +  /( y)). 

It  is  easy  to  show  that  if  the  polynomial  /  is  convex 
over  the  tctinhcdroii  I’,  then  the  implicit  patch  defined 
as  the  poition  of S[f)  inside  l’  is  convex.  However,  the 
convene  is  not  tine. 

Motivated  In  the  design  of  paianictllc'  convex  sur¬ 
faces.  icscauheis  in  CAGD  have  obtained  many  results 
on  the  convexity  conditions  of  poly  noiuials  over  trian¬ 
gles  [CPS  1].  It  is  possible  to  generalize  these  results 
to  polynomials  ovei  tetrahedra  [DMss].  thus  obtaining 
sufficient  conditions  for  implicit  patches  to  be  convex. 
Ilowevei.  the  convexity  conditions  obtained  this  way  are 
often  oveily  icstnci •'■<>.  So  in  the  following,  we  deuve 
the  convexity  conditions  of  tin  implicit  patch  directly 
fiom  the  definition  of  a  convex  implicit  patch. 

Let  p'  =  (r,'.  rj.  rj.  rj  J  be  a  point  close  to  a  point  p  = 
1  "1 .  ,  "t.  r(  |  on  I  he  sin  face  .S(  / ).  The  Taylor  expansion 

of  /.  with  higher  oidei  temis  omitted,  is 

/(p')  =  /lp>  +  ]T^ri>-, + 


Figure  3:  Weight  function*  foi  u  u ni v;i i la t o  cubic  function 
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The  first  term  on  the  right-hand  side  of  the  above 
equation,  /( p),  vanishes  because  p  is  on  $(f).  Moving 
on  to  the  second  term,  we  notice  that  this  term  is  the 
same  as  the  left-hand  side  of  the  tangent  plane  equation 
at  p, 


i=i 


So  the  definition  of  a  convex  implicit  surface  implies 
that  S(f)  is  convex  at  p  if 


since  the  barvcentiic  coordinates  (r,)  and  satisfy 
the  const  mints 

y>,  =  l  and  y  r[  =  1. 

i=t  1=1 


To  eliminate  the  const  mint  (5).  we  substitute  -  a, 
foi  i7 1  in  (-1 ).  The  iesiill  is 
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y  <t,j(p)-7,<7j  >  0  (<>) 

.  j  =  i 

for  in  bit  i  ;ti  y  (<7| .  (7j.  <7j )  with 


£ 


dTiTj 


(r/  -  r,-)(r/  -  Tj)  >  0 


(3) 


for  all  p'.  Introducing  new  variable  cr,-  =  t[  -  r, ,  we  can 
write  (3)  as 


Y'  o  ) 
dnd 


i,;=t 


d7f 

drib~]<T'°3  “  ° 


(■I) 


The  cr’s  satisfy  the  constraint 
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The  condition  (<>)  is  the  condition  foi  the  surface  S{f) 
to  be  convex  at  the  point  p. 

Applying  the  condition  (<>)  to  eveiv  point  on  an  im¬ 
plicit  patch,  we  have  the  following  theorem. 

Theorem  1  .In  implicit  / mlch  IS  (Oiiri  i  if  tin  3  x  3 
malm  .1  =  (</,_, )  is  pa nhvt  ilifinili  foi  all  po ml*  p  oil 
III <  implicit  palcli. 


i=l 


=  0 


Pioof:  Obvious  fiom  the  above  argumenls.  «!» 
(fi-neiali/.ing  the  convexity  conditions  for  bivariate 
polynomials  would  give  a  sufficient  condition  requiring 
the  matiix  .1  to  lie  positive  definite  ovei  the  entile  tetra¬ 
hedron  as  opposed  to  the  implicit  patch  The  condition 
in  Theoiein  1  is  much  less  restrictive. 


Although  Theorem  1  gives  a  condition  for  the  con¬ 
vexity  of  an  implicit  patch,  the  condition  is  hard  to  use 
because  checking  the  condition  for  the  infinitely  many 
points  on  the  implicit  patch  is  impossible.  So  in  the  rest 
of  the  section,  we  use  Theorem  1  to  derive  the  convex¬ 
ity  condition  of  an  implicit  patch  in  terms  of  its  control 
points. 

If  /  is  a  degree  k  polynomial  given  by 


/=  Y, 

m--=* 


then 


ay 

dndTj 


k(k  —  l)  bfl+e.+c ,B^  *(r), 

M=k-2 


wheie 

4 

Qm  ( o )  —  ^  ^  ( b,  I  ^  ^  ,,  m  -f 

i  j=l 

Since  the  left-hand  side  of  inequality  (11)  is  a  convex 
combination  of  Qm(cr),  the  inequality  (11)  is  \ahd  over 
the  entire  tetrahedion  enclosing  the  cubic  patch  if  and 
only  if  the  inequality  is  valid  at  the  veitices  of  the  tetra¬ 
hedion.  i.e. 


Q„ i(<7)  >0.  for  in  =  1.2.3.-1. 

So  the  cubic  patch  is  convex  if  the  inequalities  in  (8) 
holds  foi  in  =  1.2.  3.-1  with  constant 

tl,j  +  ‘V  J.,  —  h„-l.r,iT,.m  —  + 


and  A  =  (a (p))  is  a  3  x  3  symmetric  matrix  whose 
entries  are  homogeneous  degree  k  —  2  poly  nominis  in 
(r,).  From  linear  algebra,  A  is  positive  definite  if  and 
only  if 


»u  >  0, 


<Z|I<lt2 

U|2«22 


>  0.  and  |A!  >  0. 


(?) 


In  order  to  decide  whether  A  is  positive  definite  for  all 
points  p  on  an  implicit  patch,  we  have  to  determine  the 
signs  of  the  minimum  values  of  the  quantities  listed  in 
(8)  under  the  following  the  constraints, 

/(P)  =  0,  (9) 


n  +  r-i  -t- r3  +  =  1 ,  and  r,  >  0  (i  =  1,2,3, 4).  (10) 

Here  the  constraints  characterize  the  points  on  the  im¬ 
plicit  patch.  The  inequality  constraints  and  nonlinear 
constraints  in  (9)  and  (10)  make  the  problem  of  deciding 
the  convexity  of  a  general  implicit  patch  very  hard. 

Fortunately,  practical  criterions  for  the  convexity  of 
quadric  patches  and  cubic  patches  can  be  derived.  For 
quadric  patches,  notice  that 


ay 

dndTj 


=  b 


is  independent  of  p,  so  the  convexity  of  a  quadric  patch 
can  be  decided  by  evaluating  (8)  with  the  constant  ii,j  = 
6e>+c;  +  60002  —  —  bcJ+,t. 

Deciding  the  convexity  of  a  cubic  patch  is  a  little  bit 
harder.  Using  the  formula  for  Bernstein-Bezier  polyno¬ 
mials,  it  is  easy  to  verify  that 


02f 

OtiDtj 


- 6  X/  b'm r’ 


Using  this  relation,  we  can  re-  rite  (6)  as 


An  important  observation  is  that  for  each  in.  the 
alum  condition  is  exactly  the  same  as  the  convexity 
condition  foi  a  quadric  patch.  Using  the-  t<  rniiiiology  of 
(,'A(!D.  we  can  say  that  a  cubic  patch  inside  a  tetrahe¬ 
dion  is  convex  if  the  siibpolynomuiK  at  the  vertices  of 
the  teliahedion  ate  convex. 
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ABSTRACT 

Modeling  and  simulating  collections  of  physical  objects  which  are  subject  to  a  wide 
variety  of  physical  forces  and  interactions  is  exceedingly  difficult.  The  construction  of 
a  single  simulator  capable  of  dealing  with  all  possible  physical  processes  is  completely 
impractical  and.  it  seems  to  us.  wrong-headed.  Instead,  we  propose  to  build  custom 
simulators  for  single,  particular  collections  of  physical  objects  and  where  pre-specified 
physical  phenomena  are  involved.  For  such  an  approach  to  be  practical,  an  environment 
needs  to  be  provided  that  facilitates  the  quick  construction  of  these  simulators.  In  this 
paper  we  describe  the  essential  features  of  such  an  environment  and  describe  in  some 
detail  how  a  general  implementation  of  the  weighted  residual  method,  one  of  the  more 
general  classes  of  of  numerical  integration  techniques,  can  be  used. 


Keywords:  Simulation,  physical  modeling,  computational  Quid  dynamics,  symbnlii  >  .im¬ 
putation,  weighted  residual  method,  software  development  tools. 


1  Introduction 

We  are  interest?  '  in  building  software  systems  that  simulate  reality — especially 
when  several  different  physical  phenomena  are  involved  in  the  simulation.  Depend¬ 
ing  on  the  nature  of  the  objects  in  a  scene,  their  behavior  may  be  governed  by 
rigid  body  dynamics,  fluid  flow,  quantum  mechanics  or  other  families  of  laws.  The 
forces  that  act  on  these  objects  are  gravitation  and  electromagnetism  for  macro¬ 
scopic  systems,  and  weak  and  strong  interactions  for  systems  at  atomic  scales.  In 
addition,  many  observable  properties  of  physical  systems,  including  superconduc¬ 
tivity,  semiconductors  and  chemistry,  are  manifestations  of  statistical  averages  of 
detailed  lower  level  behavior.  These  macroscopic  phenomena  are  usually  simulated 
through  their  own  models — it  being  prohibitively  expensive  to  simulate  from  first 
principles. 

Besides  the  computational  costs,  the  complexity  of  dealing  with  all  physical  phe¬ 
nomena  and  mechanisms  would  make  such  a  simulator  ferociously  difficult  to  build. 
Rather  than  build  such  a  general  purpose  simulator  we  propose  a  new  program: 
Build  special  purpose  simulators  tuned  for  a  particular  configuration  of  physical 
objects  and  where  a  particular  set  of  physical  phenomena  are  involved.  Such  a 
simulator  should  be  less  complex  than  a  general  purpose  simulator,  which  must 
be  prepared  for  any  eventuality.  The  specialized  simulator  will  only  have  to  con¬ 
sider  a  known  system  of  equations  with  known  parameters.  It  will  consist  of  more 
straight-line  code  and  will  have  fewer  parameters  and  thus  should  be  easier  to  tune 
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Figure  1:  Simulation  Architecture 


for  high  performance/parallel  computer  architectures.  However,  each  new  problem 
configuration  would  require  the  creation  of  a  new  simulator. 

To  make  this  endeavor  practical,  we  are  combining  a  wide  array  of  techniques 
from  artificial  intelligence,  computer  algebra,  compiler  technology  and  code  trans¬ 
forms  to  provide  an  environment  that  vastly  simplifies  the  process  of  building  special 
purpose  simulators.  In  effect,  one  builds  a  “simulator  generator”  that  crafts  a  cus¬ 
tom  simulator  for  a  particular  configuration  of  physical  objects.  Such  a  simulator 
generator  will  generate  the  particular  set  of  differential  equations  that  model  the 
scene  and  then  will  generate  a  piece  of  code  for  the  explicit  equations  that  apply  to 
the  problem.  This  approach  has  a  number  of  advantages: 

•  More  sophisticated  mathematical  techniques  can  be  used  to  generate  the  sys¬ 
tems  of  equations  to  be  solved. 

o  Conformal  mapping  techniques  can  be  applied  to  the  non-linear  differ¬ 
ential  equations  to  simplify  and  regularize  boundary  conditions. 

o  Averaging  and  perturbation  techniques  can  be  applied  to  reduce  the 
order  of  the  equations. 

•  Numerical  techniques  specialized  to  the  equations  being  solved  can  be  used. 

•  Software  can  be  retargeted  to  different  computer  architectures  relatively  easily. 

This  new  approach  to  simulation  and  modeling  is  replete  with  new  problems  that 
need  to  be  studied  and  new  technologies  that  need  to  be  developed  and  applied.  In 
this  paper  we  sketch  a  general  framework  for  simulator  generation  and  consider  some 
of  the  components  in  detail.  It  should  be  noted  that  we  are  sketching  a  simulation 
and  analysis  framework  that  is  to  be  used  not  only  for  Newtonian  mechanics  but 
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also  for  problems  that  are  driven  by  electrodynamics,  relativistic  mechanics  and/or 
quantum  mechanics  as  well  as  aggregate  models  like  solid  state  theory,  galactic 
dynamics,  chemical  kinetics  and  fluid  dynamics.  Thus  one  should  exercise  caution 
when  extrapolating  from  experience  in  just  one  simulation  domain. 

The  process  of  performing  a  simulation  is  shown  in  Figure  1.  We  begin  with 
the  observable  scene  to  be  simulated.  An  observable  scene  is  those  properties  of 
the  system  that  ba  be  observed  and  are  independent  of  the  physics  used  to  model 
the  behavior  of  the  scene.  By  applying  the  laws  of  physics  to  the  scene  state  equa¬ 
tions  are  generated  whose  solution  describes  the  evolution  of  the  scene  with  time, 
within  certain  regions  of  validity.  The  state  equations  are  then  converted  into  code 
that  numerically  computes  the  scene’s  state  changes.  As  time  advances,  the  state 
equations  may  cease  to  be  valid  and  must  be  changed.  Similarly,  the  geometric  or 
topological  characteristics  of  the  scene  itself  may  change.  These  effects  are  indicated 
by  the  shaded  “feedback”  arrows  in  Figure  1. 

The  state  of  a  physical  system  is  determined  by  the  values  of  a  set  of  state 
variables,  which  may  include  a  subset  of  the  observable  parameters  of  the  objects 
in  a  scene.  The  result  of  applying  the  physical  laws  to  a  scene  are  a  set  of  stale 
equations  that  constrain  the  state  variables  over  time.  For  instance,  for  clocked 
boolean  logic  circuits,  there  is  a  finite  set  of  state  variables,  each  of  which  ranges 
over  {true,  false},  time  is  modeled  by  a  sequence  of  discrete  events  occurring  on 
clock  edges  and  the  state  equations  are  boolean  equations.  For  rigid  body  dynamics, 
there  is  a  discrete  set  of  state  variables  that  have  continuous  values,  time  is  modeled 
as  a  sequence  of  irregularly  spaced  events  and  the  state  equations  are  ordinary 
differential  equations.  In  fluid  dynamics,  there  is  a  continuous  number  of  of  state 
variables,  one  for  each  point  in  the  (continuous)  fluid,  and  their  values  range  over 
a  continuous  vector  space.  The  state  equations  are  partial  differential  equations. 

Once  the  state  equations  of  a  scene  have  been  genera' ed  (the  middle  box  of 
Figure  1),  general  mathematical  techniques  can  be  used  to  convert  them  into  a 
form  where  numerical  information  about  their  state  variables  can  be  determined. 
Examples  include  conversions  of  ordinary  differential  equations  into  finite  differ¬ 
ence  formulas  by  Runge-Kutta  methods,  or  the  conversion  of  partial  differential 
equations  into  systems  of  linear  equations  by  finite  element  methods.  We  call  the 
process  of  converting  a  system  of  equations  into  an  effective  computational  form  a 
discretization. 

These  computation  structures  can  then  be  converted  to  actual  programs  (or 
codes)  that  numerically  simulate  the  scene.  If  something  is  known  about  the  archi¬ 
tecture  of  the  computer  that  will  run  the  program  then  especially  fast  codes  can 
be  generated  by  symbolic  elimination  of  variables,  unrolling  of  loops  or  duplication 
of  code.  Each  of  these  options  may  be  appropriate  because  of  cache  sizes,  vec¬ 
tor  processing  structures  or  interprocessor  communication  costs.  Other  .echniques 
of  compiler  theory  are  also  appropriate  and  should  be  carefully  considered  at  this 
point.  More  radical  transformations  like  changing  the  order  of  the  discretization 
or  the  discretization  method  itself  may  also  be  appropriate.  This  entire  process  of 
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converting  state  equations  into  computational  structures  and  then  into  executable 
code  is  indicated  in  the  right  half  of  Figure  l. 

This  paper  discusses  each  of  these  steps  in  the  simulation  process.  In  Section  2  we 
discuss  one  approach  to  representing  scenes,  their  components  and  the  underlying 
physics.  Once  the  state  equations  have  been  generated,  they  can  be  directly  solved 
numerically,  yielding  the  trajectory  of  the  scene  from  a  given  set  of  initial  conditions. 
The  approach  we  are  pursuing  is  discussed  in  some  detail  in  Section  4. 

However,  occasionally  some  property  of  the  trajectory  is  of  interest — not  the  tra¬ 
jectory  itself.  We  argue  that  by  using  symbolic  techniques,  we  can  often  transform 
the  differential  equations  that  describe  the  system  into  other  equations  whose  solu¬ 
tions  more  precisely  answer  the  questions  being  asked.  Solving  these  transformed 
equations  is  often  substantially  easier  than  solving  the  original  system.  However, 
substantial  (non-numeric  computation)  is  required  to  produce  the  transformed  equa¬ 
tions.  In  Section  3  we  illustrate  how  averaging  techniques  can  be  applied  to  reduce 
the  dimension  of  the  problem  being  solved  and  more  directly  answer  the  questions 
of  interest.  This  technique  is  classical,  but  we  fell  is  representative  of  the  type  of  re¬ 
duction  that  will  be  valuable  in  the  future  and  is  possible  by  the  general  framework 
being  proposed. 

In  the  domain  in  which  we  are  working  (fluid  dynamics),  the  state  equations 
are  partial  differential  equations.  A  wide  variety  of  different  methods  are  available 
for  their  numerical  solution.  Many  of  these  methods  can  be  subsumed  within  the 
general  mathematical  framework  of  weighted  residual  methods.  Because  we  have 
access  to  the  state  equations  in  symbolic  form,  we  can  directly  apply  the  weighted 
residual  methodology  to  the  differential  equations  of  the  problem  to  produce  a 
computational  structure  based  on  a  wide  variety  of  different  techniques  including 
finite  element,  spectral  and  collocation  methods.  This  approach  is  discussed  in 
Section  4.  In  Section  4.1  we  describe  the  general  principles  behind  the  weighted 
residual  method.  In  Section  4.2  we  use  the  weighted  residual  method  to  produce 
a  spectral  method  computational  structure  for  a  problem  in  fluid  dynamics.  This 
particular  example  illustrates  the  complexity  of  the  codes  generated  in  the  study  of 
turbulent  fluid  dynamics. 

In  Section  4.3  we  give  another  illustration  of  the  weighted  residual  method  in 
fluid  dynamics,  but  this  time  the  result  of  the  discretization  process  is  not  a  system 
of  linear  equations,  but  rather  a  system  of  ordinary  differential  equations.  This  is 
another  example  of  where  symbolic  techniques  can  be  used  to  convert  a  numerical 
problem  into  one  that  more  directly  provides  the  desired  answers. 

2  Scenes  and  Laws  of  Physics 

This  section  makes  more  precise  what  we  mean  by  scenes  and  physical  laws.  Sec¬ 
tion  2.1  discusses  scenes  while  Section  2.2  discusses  the  components  of  a  physics 
and  some  their  functions. 
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2.1  Scenes 


When  describing  a  physical  system  that  is  to  be  simulated,  we  distinguish  the 
observable  properties  and  characteristics  of  the  system  from  those  properties  and 
characteristics  that  are  required  by  a  particular  physical  model.  The  former  are 
components  of  the  observable  scene,  while  the  later  belong  to  the  physical  laws 
and  models  that  are  to  be  applied  to  the  scene.  For  instance,  the  charge,  mass 
and  position  of  an  electron  are  components  of  the  scene,  but  the  electric  field  is  a 
characteristic  of  a  physical  model  that  might  be  used  to  determine  the  efTect  of  the 
electron  on  other  charged  particles.  Fields  in  physics  are  not  themselves  directly 
observable.  It  is  only  through  their  effect  on  other  objects  that  their  existence  can 
be  ascertained.  The  effect  of  one  object  on  another  is  the  purpose  of  a  physics 
and  thus  fields  are  artifices  used  to  facilitate  the  physics  itself.  (Recall  that  general 
relativity  replaces  a  gravitational  field  by  bending  of  the  fabric  of  space  itself.  For 
small  masses  these  disparate  mechanisms  give  the  same  predications.) 

Similarly,  the  “ether”  of  nineteenth  century  physics  belongs  to  a  set  of  physical 
laws,  and  is  not  intrinsic  to  the  scene.  Ether  is  posited  by  nineteenth  century 
physics  and  is  not,  itself,  observable.  A  more  modern  example  is  the  wave  function 
of  quantum  mechanics.  It  cannot  be  observed  in  the  scene  but  is  essential  for  a 
particular  set  of  physical  laws.  In  all  of  these  cases  the  physics  used  to  analyze  the 
system  imposes  additional  parameters  ( e.g .,  wave  functions)  or  objects  ( e.g .,  fields 
or  ether)  as  an  aid  in  specifying  the  physics  itself. 

A  scene  consists  of  a  number  of  objects  (rods,  resistors,  fluids,  etc.)  and  connec¬ 
tions  (hinges,  electrical  nodes,  etc.)  between  them.  The  constraints  constrain  the 
behavior  of  two  or  more  objects  in  some  fashion.  For  instance,  a  hinge  between  two 
rods  requires  that  the  rods  remain  connected,  while  an  electrical  node  connecting 
the  pins  of  two  resistors  ensures  that  the  two  pins  always  have  the  same  potential. 

In  addition,  the  objects  possess  a  number  of  “observable”  properties,  e.g.,  the 
position  and  momentum  of  a  particle.  These  properties  are  those  aspects  of  the 
state  of  the  object  that  may  be  observed  in  the  scene,  and  thus  are  independent 
of  the  physics  used  to  model  the  behavior  of  the  scene.  The  observable  properties 
may  be  redundant  and  related  by  some  equations.  For  instance,  the  observable 
properties  of  a  particle  include  the  particle’s  mass  (m),  position  (r).  velocity  (v). 
momentum  (p)  and  kinetic  energy  (T),  where 


dr 


p  =  mv, 

_  mv  -v  _  p  p 
2  2m  ' 

For  some  models  of  physics,  like  Newtonian  mechanics,  the  observable  parameters 
actually  correspond  to  state  of  the  physics.  That  is,  the  observable  position  and 
momentum  are  actually  the  position  and  momentum  of  the  object  in  the  physics. 
In  other  models,  e.g.  quantum  mechanics,  the  observables  are  derived  from  the 
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Figure  2:  Simple  Pendulum 


their  correspondents.  That  is  the  quantum  mechanical  position  and  momentum  of 
a  particle  are  not  interchangeable  with  the  observable  position  and  momentum  of 
the  particle. 

2.2  Physics 

The  properties  of  an  observable  scene  are  not  necessarily  appropriate  for  simula¬ 
tion.  Instead,  the  physical  laws  translate  the  scene  into  one  where  the  new  scene's 
objects  are  described  using  state  variables.  For  instance,  a  two  dimensional  scene 
that  consists  of  a  heavy  bob  at  the  end  of  a  rigid,  massless  rod  whose  other  end  is 
hinged  (i.e..  a  two  dimensional  pendulum)  might  have  constitutive  parameters  of 
the  length  of  the  rod  (()  and  the  mass  of  the  bob  (m) — see  Figure  2.  The  observable 
parameters  in  the  scene  might  be  the  position  of  the  bob  ({r.  y)  6  ?.2).  However, 
when  formulating  a  simulation,  one  would  probably  use  the  deflection  of  the  rod 
from  vertical  ( 0  €  [— x,  x))  as  the  state  variable  of  the  system.  The  position  of  the 
bob  can  be  derived  from  6  by 

(x,y)  =  (c-  -r  £sin0,c7  -f  (cos 9). 

Each  set  of  physical  laws  acts  similarly.  It  must  construct  from  an  observable 
scene  an  interpretation  scene  that  consists  of  objects,  state  variables  that  are  appro¬ 
priate  to  the  physical  model  and  the  manifold  structure  in  which  the  state  variables 
lie.  A  correspondence  also  needs  to  be  provided  between  state  variables  in  the  in¬ 
terpretation  scene  and  quantities  in  the  observable  scene.  The  combination  of  the 
state  variables,  their  manifold,  and  the  correspondence  we  call  an  interpretation 
scene  or  just  an  interpretation.  Examples  of  interpretations  are  the  generalized 
coordinates  of  Hamiltonian  mechanics  (which  were  used  in  the  pendulum  example) 
and  the  wave  functions  of  quantum  mechanics. 
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A  description  of  a  physics  consists  of  (1)  the  domains  of  validity  of  the  physics, 
(2)  a  method  to  generate  an  interpretation  scene  from  an  observable  scene  and  (3) 
how  fields  and  energies  are  to  be  derived  from  the  resulting  interpretation  scene.  In 
addition  we  must  have  a  formulation  of  mechanics  that  allows  us  to  combine  the 
fields  and  energies  produced  by  the  different  sets  of  physics  to  generate  the  state 
equations.  Examples  of  formulations  are  Lagrangian  and  Hamiltonian  mechanics, 
both  of  which  can  be  applied  outside  the  domain  of  Newtonian  mechanics. 

The  physical  laws  that  apply  to  a  scene  are  kept  separate  from  the  scene  and 
should  be  expressed  independently  of  their  application  to  a  particular  scene.  There 
should  be  one  (or  more)  descriptions  of  rigid  body  dynamics  and  one  (or  more) 
descriptions  of  electrodynamics.  These  descriptions  include  the  “laws  of  physics" 
( e.g .,  F  =  ma  for  rigid  body  dynamics,  or  Maxwell’s  equations),  specifications  of 
when  the  particular  laws  are  applicable  and  procedural  specifications  of  how  to 
apply  the  laws  to  a  particular  scene. 

We  call  a  set  of  physical  laws  a  physics.  Each  physics  has  a  limited  range  of  ap¬ 
plicability  (until  the  Grand  Unified  Theory  is  discovered).  Among  the  components 
of  a  physics  is  a  specification  of  how  forces  and  energies  of  objects  in  a  scene  can  be 
computed.  There  are  multiple  physics’s,  some  of  which  are  compatible  with  each 
other  over  certain  ranges  of  state  variables  and  some  of  which  are  incompatible. 
For  instance,  Newtonian  ’chanics  and  classical  electrodynamics  are  compatible 
for  small  masses  and  slow.,  noving  macroscopic  particles.  Electrodynamics  merely 
introduces  a  new  force,  which  is  characterized  by  Coulomb’s  law.  Quantum  me¬ 
chanics  and  general  relativity  seem  incompatible. 

This  approach  ensures  that  different  physical  considerations  are  dealt  with  sep¬ 
arately.  For  instance,  one  should  be  able  to  simulate  an  electric  motor  by  applying 
both  rigid  body  dynamics  and  electromagnetics  to  a  scene  that  consists  of  the  rotor 
and  stator  of  the  motor,  with  the  appropriate  constitutive  properties.  We  believe 
that  greater  modularity  will  result  from  this  approach,  althougn  it  places  a  premium 
on  the  symbolic  techniques. 

3  Harmonic  Balance 

When  setting  up  a  system  of  differential  equations  that  models  some  physical  sit¬ 
uation,  it  is  often  easier  to  generate  the  equations  in  terms  of  state  variables  that 
are  different  from  the  ones  that  the  user  is  really  interested  in.  For  instance,  for 
a  mechanical  system  it  may  be  easiest  to  generate  equations  in  terms  of  cartesian 
coordinates  while  the  interesting  behavior  might  best  be  expressed  in  terms  of  radial 
coordinates,  or  angular  momenta  or  even  averaged  angular  momenta.  Each  of  these 
conversions  can  be  performed  after  the  numerical  solutions  are  generated,  however, 
using  symbolic  techniques  to  perform  this  conversion  before  the  integration  process 
makes  generating  an  accurate  solution  easier. 

To  illustrate  this  approach  we  will  use  a  more  sophisticated  type  of  coordinate 
change  that  '.Iso  facilitates  an  averaging  technique.  Thus  we  will  ultimately  generate 
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differential  equ.  ons  for  the  average  values  of  the  variables  of  interest. 

A  large  variety  of  simple  oscillatory  type  systems  can  be  modeled  by  differential 
equations  of  the  form 

i  =  y, 

,,  ^  (1) 

y  =  -x+eh(x,y). 

For  h(x,  y)  =  (1  ~x2)y  we  have  the  van  der  Pol  equation  [14],  for  h(x ,  y)  =  (1  - y2)y 
the  Rayleigh  equation  [12],  etc.  When  c  =  0  (1)  reduces  to  a  simple  harmonic 
oscillator,  whose  solution  is: 

x(t)  =  r0  cos (t  +  <f>o),  y{t)  =  i(ar)  =  -r0  sin(f  +  <£0),  (2) 

where  ro  and  <j>o  are  constants  set  by  the  initial  conditions.  In  the  x-y  plane  (the 
phase  plane),  the  solutions  are  circles  centered  at  the  origin.  The  term  ch(x,y)  of 
(1)  acts  as  a  perturbing  non-linear  damping  factor  on  the  solution  to  the  harmonic 
oscillator.  An  example  of  the  behavior  of  this  damping  factor  can  be  seen  from  the 
van  der  Pol  equation  where  h(x,y)  =  (1  -  x2)y: 

x  =  y, 

(3) 

y  =  -x  +  c(l  -  1 2)y. 

The  phase  plot  of  the  van  der  Pol  equation,  for  c  =  0.6  and  various  initial  conditions, 
is  shown  in  Figure  3. 

In  the  phase  plane,  (3)  has  a  stable  limit  cycle  of  radius  approximately  2.  If  the 
initia'  conditions  of  the  system  are  outside  the  limit  cycle,  the  system  will  cycle 
inwards  around  the  limit  cycle  continuously  getting  closer.  If  the  starting  point  is 
inside  the  the  limit  cycle  the  system  will  oscillate  outwards  towards  the  limit  cycle. 
From  a  physical  point  of  view  we  mhjht  have  two  basic  questions: 

•  What  is  the  average  amplitude  of  the  Pratt  cycle? 

•  How  quickly  does  the  system  converge  to  the  'ira.it  cycle? 

We  can  study  the  behavior  of  the  non-linear  oscillator  by  ass.,n-ing  the  ,?olu'.:'>i 
is  of  the  form  (2)  but  allow  the  constants  to  be  time  varying  functions,  :  . 

x  =  r(f)cos(t  +  <f>{t)), 

y  =  r(t)sin(t  +  <j>(t)). 

Substituting  these  expressions  into  (1)  gives  the  following  system  of  equations 
rcos(t  +  <j>)  -  sin(<  +  <£)(1  +  j>)  =  rsin(<  +  <t>), 
rsin(t  +  <j>)  +  cos(t  +  <p)(l  4-  <j>)  =  eh(r  cos(<  +  <p),  rsin(Z  +  4>)). 

When  solved  for  r  and  0,  which  must  be  done  symbolically,  we  have 
r  =  eh(rcos(t  +  <f>),rs\n(t  +  <£))sin(<  +  <j>) 

j>  =  -^A(rcos(Z  +  <f>),  rsin(<  +  <f>))cos(t  +  <f>). 
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Figure  3:  van  der  Pol  Oscillator 


The  solution  to  this  system  of  equations  gives  the  amplitude  of  the  system,  which 
is  closer  to  what  we  are  looking  for.  In  a  physical  system  we  probably  don’t  care 
about  the  phase  information.  We  are  more  interested  in  the  asymptotic  behavior  of 
the  system.  This  can  be  obtained  by  averaging  these  equations  over  one  oscillation, 
i.e.  t  to  t  +  2 7r: 

2x 

h(rcos(t  +  <j>),  rsin(t  +  0))sin(<  +  <j>)dt, 

r2* 

—  /  h(r  cos(t  +  4>),  rsin(<  +  $))  cos(t  +  4>)dt. 

dt  2itr  Jo 

In  the  r-<f>  coordinate  system,  the  van  der  Pol  equation  becomes 


d(r)  _  c  f 

dt  2t  J0 


r  =  cr(l  —  r2  cos2(t  +  <^))sin2(t  +  <j>) 

4>  =  e(I  —  r2  cos2(<  +  <j>))  sin(i  +  <j>)  cos  (t  +  <f>) 


When,  averaged,  the  equation  for  r  becomes 


d(r) 

dt 


[2* 

—  /  (r2cos2{/  <f>)  -  l)rsin2(t  +  <f>)dt  =  -e 

2jt  Jo 


(4) 


(5) 


The  solution  of  this  equation  is  precisely  the  evolution  of  the  “average”  amplitude 
of  the  oscillation  without  any  additional  information.  Notice  that  we  have  been  able 


9 


0  2  4  6  8  10  12  14 


Figure  4:  Amplitude  of  van  der  Pol  Oscillator:  raw  and  averaged 

to  reduce  the  order  of  the  equation  by  one  by  averaging  out  the  phase  information. 
In  Figure  4  we  have  shown  the  evolution  of  the  a  solution  of  (4)  for  two  starting 
points,  one  inside  and  one  outside  the  limit  cycle,  using  a  solid  line.  The  dotted 
lines  indicate  solutions  of  the  averaged  equation  (5)  from  the  same  two  starting 
points. 

On  the  stable  limit  cycle  of  the  system,  (r)  will  vanish,  so  by  solving 


we  see  that  the  average  radius  of  the  limit  cycle  is  2,  which  is  independent  of  the 
initial  conditions.  This  can  also  be  observed  from  Figure  4. 

The  rate  at  which  a  solution  approaches  the  limit  cycle  can  be  determined  by 
solving  (5): 


This  type  of  perturbation  analysis  has  been  used  by  in  celestial  mechanics  since 
the  time  of  Laplace  and  Lagrange.  The  particular  problem  we  consider  here,  the 
behavior  of  solutions  of  equations  of  the  form  (1)  was  discussed  in  some  detail  by 
Poincare  [11].  More  recently  Krylov  and  Bogoliubov  [8]  have  demonstrated  the  use 
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of  averaging  techniques  in  a  wide  variety  of  problems.  For  a  more  modern  treatment, 
one  might  look  at  [13]. 

4  Weighted  Residual  Methods 

For  a  large  number  of  physical  simulation  problems  the  state  equations  are  partial 
differential  equations.  Though  the  number  of  techniques  for  solving  these  systems 
can  be  bewildering  in  their  number,  the  most  important  techniques  can  be  divided 
into  two  major  classes:  finite  difference  algorithms  and  weighted  residual  meth¬ 
ods.  We  have  decided  to  focus  on  weighted  residual  methods  because  most  of  the 
techniques  of  interest  in  our  application  area  are  of  the  weighted  residual  type. 

There  are  a  vast  number  of  implementations  of  numerical  algorithms  based  on 
particular  weighted  residual  methods,  most  often  for  particular  partial  differential 
equations,  but  to  our  knowledge  there  have  been  no  previous  attempts  to  build 
a  system  that  generates  a  numerical  solver  for  a  wide  class  of  weighted  residual 
methods. 

We  describe  the  basic  principles  behind  the  weighted  residual  method  in  Sec¬ 
tion  4.1.  In  Section  4.2  we  give  a  brief  illustration  of  how  the  weighted  residual 
method  is  used  to  generate  particular  numerical  codes  in  fluid  mechanics.  Finally, 
in  Section  4.3  we  use  the  weighted  residual  method,  along  with  a  number  of  other 
ideas,  to  reduce  some  questions  about  fluid  flow  to  questions  about  a  system  of 
ordinary  differential  equations. 

4.1  General  Approach 

Let 

Lu  =  f  (6) 

be  a  partial  differential  equation,  where  L  is  a  partial  differential  operator  and  u  is 
a  function  of  {xi,..., xm}.  The  weighted  residual  method  assumes  there  exists  a 
(possibly  infinite)  set  of  trial  functions  {&}  such  that,  for  some  choice  a,-, 

«=  ai<pi  (7) 

0  <i<N 

is  an  approximation  to  u,  the  solution  of  (6).  The  <£,-  are  function  of  some  subset 
of  {xi, . . . ,  xm}  while  the  o;  are  functions  of  the  remaining  variables.  Substituting 
(7)  into  (6)  we  have  residual  error 


Re(u)  =  L 


-/• 


The  goal  of  a  weighted  residual  method  is  to  choose  the  a,  in  a  fashion  that  mini¬ 
mizes  Re(u)  in  some  global  sen'  . 
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A  set  of  equations  for  the  a,-  can  be  deduced  by  choosing  a  set  of  weighting 
functions,  w}  and  requiring  the  inner  product  of  Re(u)  with  the  weights  to  vanish, 
i.e. 

J  WjR.E{ii)dV  =  0  (8) 

If  the  4>i  are  functions  of  all  of  the  variables  then  the  resulting  equa¬ 

tions  are  algebraic  in  the  ar  When  I  is  a  linear  differential  operator,  the  resulting 
equations  are  linear.  Applying  I  to  the  components  of  the  expansion  and  rewriting 
(8)  we  have 

Oj  I  Wj  L<f>idV  =  f  Wjf  dV 

o  <i<N 

For  certain  operators  L  and  known  <p,  and  w}  the  integrals  above  can  be  tabulated. 
Thus  the  bulk  of  the  symbolic  computation  inherent  in  the  reduction  of  (8)  to 
systems  of  equations  in  the  a,  can  be  performed  a  priori.  However,  if  the  6,  and 
are  supplied  by  the  user  and,  especially,  if  L  is  non-linear  then  symbolic  computation 
is  unavoidable  in  the  application  of  the  weighted  residual  method. 

If  the  cpi  involve  a  subset  of  {xi,...,xm}  then  (8)  will  be  a  system  of  ordinary 
differential  equations.  Often  the  <j>,  are  functions  only  of  the  spatial  variables  and 
the  a}  are  functions  of  time.  This  is  the  situation  in  the  two  examples  considered 
here.  In  Section  4.2  the  partial  differential  equation  is  first  discretized  in  time 
and  then  the  weighted  residual  method  is  applied,  producing  a  system  of  linear 
equations  that  need  to  be  solved.  In  Section  4.3,  the  weighted  residual  method  is 
applied  directly  to  the  spatial  variables  resulting  in  a  system  of  ordinary  differential 
equations  for  the  a,-. 

A  wide  variety  of  different  integration  schemes  fall  into  this  general  framework. 
If  the  Wj  are  chosen  to  be  the  same  as  the  $,•  we  get  a  Galerkin  projection.  This  is 
especially  convenient  if  the  <f>,  are  orthogonal  and  eigenfunctions  of  L.  The  system 
of  linear  equations  are  then  diagonal.  The  resulting  technique  is  called  a  spectral 
method.  The  most  common  spectral  method  chooses  4>t  =  e'^'x. 

The  finite  element  method  discretizes  the  computational  domain  fl  into  a  number 
of  elements,  It  then  chooses  the  <£,•  to  be  continuous  functions  that  are 

zero  everywhere  except  within  fl,.  A  Galerkin  projection  then  gives  the  equations 
for  the  a,-. 

In  general,  determining  the  linear  equations  or  ordinary  differential  equations 
that  need  to  be  solved  from  (8)  is  a  rather  painful  process  that  must  be  performed 
by  hand.  By  taking  advantage  of  methods  from  symbolic  computation  we  can 
largely  automate  this  process. 

4.2  Numerical  Example 

In  this  section  we  illustrate  how  the  weighted  residual  method  is  used  to  produce 
a  numerical  code  for  a  problem  that  arises  in  the  study  of  turbulent  channel  flow. 
This  example  illustrates  the  complexities  that  arise  in  practical  applications  of  the 
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weighted  residual  method.  A  simplification  of  one  of  the  equations  that  arose  in 
Kim,  Moin  and  Moser’s  study  of  turbulence  in  a  channel  flows  [6]  is 


dg(x,y,L) 

dt 


=  Kg)  + 


(9) 


where  x  and  y  are  the  spatial  dimensions  of  the  the  problem  (only  two  are  needed 
for  this  illustration),  h  is  a  known  non-linear  function  of  g  and  other  functions  that 
occur  in  the  problem.  In  practice  it  can  be  a  fairly  complex  expression  and  more 
than  one  partial  differential  equation  be  involved. 

In  solving  this  problem,  discretization  must  occur  in  three  different  dimensions — 
time  and  the  two  spatial  dimensions.  Three  different  schemes  will  be  used:  An 
implicit  finite  difference  scheme  for  time  (t),  a  spectral  method  for  the  x  dimension 
and  a  Galerkin  type  method  using  Chebyshev  polynomials  for  the  y  dimension. 

These  three  schemes  are  used  in  three  successive  steps.  First,  time  is  discretized 
and  the  value  of  g(x,  y,  t)  at  the  nth  time  step  is  denoted  by  gn(x,  y)  =  g(x,  y,  n  At). 
Second,  gn(x,y)  is  discretized  in  the  x  dimension  using  Fourier  expansion  with 
coefficients  y£(y).  Finally,  y£(y)  is  discretized  using  Chebyshev  polynomials  in  the 
y  dimension  with  coefficient  y£ That  is, 


<?"(*,  y)=  E  9nt(y)e2*'kl,L', 

0 <k<Nt 

=  E  E  rkMy)^ikt'Lx- 

0<k<Nr0<j<Nv 


At  this  point  the  coefficients  are  numbers,  and  if  done  properly  they  are  solutions  of 
linear  equations.  Once  these  linear  equations  have  been  solved  we  can  reconstruct 
g(x,y,t)  by  summing  the  series. 

Each  of  these  transformations  can  be  automated  using  symbolic  techniques.  In 
practice,  their  application  is  not  completely  straightforward.  The  following  para¬ 
graphs  illustrate  this  with  some  comments  on  the  implementation  of  these  tech¬ 
niques  using  symbolic  computation. 

The  first  step  is  to  perform  the  time-wise  discretization.  We  denote  by  gn  = 
gn(x,y)  the  function  that  corresponds  to  g  at  the  nth  time  step.  The  most  straight¬ 
forward  discretization  would  be  the  explicit  formula 


=  W)  +  i£V’, 


n 


But  this  is  known  to  be  relatively  unstable. 

Figure  5  gives  a  number  of  different  discretization  techniques  that  can  be  used 
for  ordinary  differential  equations.  For  this  particular  problem  no  one  of  them 
is  completely  satisfactory.  For  instance,  the  explicit  methods  (Euler  or  Adams- 
Bashforth  (2  4])  are  not  sufficiently  stable  when  applied  to  the  entire  equation. 
The  implicit  methods  require  the  solution  of  a  non-linear  equation  at  each  time 
step  (because  of  the  nonlinearity  of  h )  and  are  thus  too  costly. 
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j.n+1  _  xn 


At 


'  /(*") 
f(xn+l) 

f(xn+l)+f(xn) 

*  2 

±[23f(xn)-l6f(xn-l)  +  Sf(tn-i)] 


Explicit  Euler 
Implicit  Euler 
Crank-Nicolson 
2nd  Adams- Bashforth 
3rd  Adams- Bashforth 


Figure  5:  Discretization  techniques  for  i(t )  =  /(x) 


The  solution  is  to  use  an  explicit  scheme  on  the  nonlinear  terms  and  an  implicit 
scheme  on  the  linear  terms.  Using  the  second  order  Adams- Bashforth  formula  for 
the  linear  terms  and  the  Crank-Nicolson  formula  [3]  for  the  linear  terms  yields 

=  £  (3%n)  -  %-1))  +  ok (vV+1 + vV)  • 

In  a  symbolic  manipulation  system  this  process  is  quite  simple.  The  differential 
equation  is  first  converted  to  a  sum  of  terms  form.  Each  term  is  then  examined  to 
see  if  it  is  linear  in  g.  If  so,  an  implicit  formula  is  applied  to  each  term,  otherwise 
an  explicit  one.  The  results  of  these  replacements  are  then  added  together  and 
simplified. 

The  terms  that  involve  gn+1  can  be  isolated  on  the  left  hand  side  to  give 

sn+1 "  ^vV+1  =  T (3%n)  “  /,(?n"1))  +?n  +  (10) 

Again  the  symbolic  processing  involved  is  straightforward,  each  term  is  examined 
and  placed  on  one  side  of  the  equation  or  the  other  based  on  its  dependence  on 
gn+l. 

At  a  given  time  step,  each  of  the  terms  on  the  right  hand  side  of  (10)  is  known 
and  can  be  computed  directly.  The  next  step  is  to  compute  the  Fourier  transform 
of  this  equation,  eliminating  the  functional  dependence  on  x. 

9n(x> y)  =  £  gnk(yy*ikx/Lz- 

0  <k<Nz 


Thus  the  kth  mode  the  Fourier  transform  of  the  left  hand  side  of  (10)  is 


= 1 


14 


So  for  each  k  we  need  to  solve  the  equation: 


This  equation  is  finally  discretized  in  y  by  expanding  y£(y)  in  terms  of  Chebyshev 
polynomials: 

fi(y)=  £  KjTiiy), 

0  <}<NV 

where  y£ ;  are  numbers.  The  Chebyshev  expansion  of  the  left  hand  side  of  (11)  is 


The  last  term  in  this  sum  causes  some  problems  because  it  is  not  expressed 
as  a  sum  of  Chebyshev  polynomials,  but  as  a  sum  of  their  derivatives.  However, 
derivatives  of  Chebyshev  polynomials  can  be  expressed  as  a  sum  of  Chebyshev 
polynomials  by  repeated  application  of  the  formula 

^(*)  to  .  ??-,(«)  _1T  (x) 

(n+l)(n  +  2)  (n2  —  1)  (n  —  l)(u  —  2)  h 

or  by  solving  the  tridiagonal  system  it  implies.  At  this  point,  we  have  converted 
the  problem  of  advancing  time  in  (9)  to  solving  systems  of  linear  equations  and 
computing  Fourier  and  Chebyshev  transforms. 

For  other  basis  and  weight  functions,  and  for  other  differential  equations,  very 
similar  approaches  are  used.  Simple  symbolic  methods  (arithmetic  operations  and 
some  simplification)  are  used  to  reduce  the  projection  process  to  a  sequence  of 
integral.  In  the  case  discussed  here,  all  of  the  integrals  could  be  performed  by  table 
lookup.  In  the  next  section  the  integrals  will  have  to  be  performed  numerically. 

4.3  Proper  Orthogonal  Decomposition 

By  discretizing  the  spatial  dimensions  but  not  the  time  dimension,  we  can  reduce  the 
Navier-Stokes  equations  to  a  system  of  ordinary  differential  equations.  If  the  proper 
basis  functions  are  chosen  and  sufficient  terms  are  used  the  dynamical  behavior  of 
the  ODE’s  should  closely  approximate  that  of  the  Navier-Stokes  equations. 
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Figure  6:  Coordinate  System  for  a  Channel 


Lumely  [9]  has  suggested  using  this  approach  to  study  the  behavior  of  the  turbu- 
len«  boundary  layer  of  a  fluid  moving  over  a  flat  plate.  Within  the  boundary  layer 
bursts  can  be  observed  that  are  spatially  and  temporally  somewhat  periodic  [7].  It 
would  be  interesting  to  know  if  these  periodic  phenomena  manifests  themselves  in 
the  ordinary  differential  equations  where  more  powerful  mathematical  techniques 
can  be  used  to  analyze  their  behavior. 

This  reduction  and  detailed  study  of  the  resulting  dynamical  systems  was  origi¬ 
nally  performed  by  John  Lumley,  Philip  Holmes  and  their  students  (I).  As  we  shall 
see  the  ordinary  differential  equations  that  result  are  extremely  complex  and  are 
best  generated  by  symbolic  techniques. 

Fluid  flow  is  governed  by  the  Navier-Stokes  equations.  In  the  absence  of  external 
forces  the  dimensionless  form  of  these  equations  is 

^ +  (v  •  V)v  = -Vtt  + ■i-V2v,  (12) 

V  •  v  =  0  (13) 

where  v  denotes  the  velocity  field  of  the  fluid  and  Re  is  the  Reynolds  number  of 
the  fluid.  Flows  with  small  Reynolds  numbers  tend  to  be  rather  steady,  while  flows 
with  Reynolds  numbers  greater  than  2300  are  generally  turbulent.  In  order  to  write 
the  equations  in  a  dimensionless  form,  characteristic  lengths  need  to  be  defined  in 
each  of  the  three  dimensions.  We  denote  these  different  characteristic  lengths  by 
£i,  L2  and  £3. 

If  f(xi,x 2,13)  is  function  of  position  in  the  channel,  we  will  denote  by  (/)  its 
spatial  average  in  a  plane  parallel  to  the  walls  of  the  channel: 

(f(xi,x2,z3))  =  -  J  f(x  1 ,  x2,  x3)  dx \dx3, 


16 


where  Ly  and  I3  are  characteristic  dimensions  in  the  xy  and  £3  directions  respec¬ 
tively.  Within  this  plane  the  turbulent  flow  is  relatively  homogeneous.  Variations 
occur  in  the  orthogonal  direction.  Thus  (/(ij,  x2,  X3))  is  a  function  of  only  the 
distance  from  the  wall,  x2. 

The  streak  structure  that  we  are  interested  in  is  not  a  function  of  the  mean 
velocity  of  the  fluid,  only  its  fluctuations.  Denoting 

(v)  =  (f/(x2),  0,0)  =  U, 
we  can  determine  U(x 2)  exactly: 

U{x2)  =  ReJ  (uyu2)dx'2  +  Reu£  ^x2  -  ^  .  (14) 

where  ut  is  the  dimensionless  wall  shear  velocity  and  H  is  the  half  height  of  the 
channel.  Applying  this  to  (12)  we  get  the  Reynolds  averaged  Navier-Stokes  equa¬ 
tions: 


(15) 

These  are  the  equations  to  which  we  will  apply  the  weighted  residual  method. 

Notice  that  while  the  Navier-Stokes  equations  are  quadratic  these  equations  for 
the  velocity  fluctuation  are  cubic  due  to  the  quadratic  behavior  of  the  mean  velocity 
in  (14). 


4.3.1  Eigenfunction  Projections 

Because  the  flow  is  homogeneous  in  the  plane  parallel  to  the  wall,  we  can  use  a 
Fourier  expansions  in  the  in  xy  and  £3  directions  (parallel  to  the  wall).  We  assume 
we  are  given  a  set  of  basis  functions  in  the  inhomogeneous  direction.  These  basis 
will  not  be  known  numerically,  but  rather  will  be  provided  numerically.  These 
basis  functions  are  called  the  “empirical  eigenfunctions.”  Expanding  the  velocity 
fluctuation  u  just  along  Xy  and  £3  dimensions  we  have 

u(xi,x2,x3,t) 

fcj  =  -CO 

Each  of  the  ii(x2,f,ky,k3)  can  be  expanded  in  series  based  on  the  empirical  eigen¬ 
functions: 

CO 

u{x2)t\ky,k3)  = 

n  =  1 


VLIL3 


f;  u(£2,t;*l,*3)e2T,<^r,+£r3). 
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Combining  these  two  expansions  gives  the  following  representation  of  the  \elocity 
fluctuation  field: 


.  CO  o° 


(16) 


Jt3=-oo 


Notice  that  (16)  is  actually  three  equations,  one  for  each  component  of  u.  We 
denote  the  components  of  <4"*,  by 


4> 


(") 


,4> 


(") 


.(» 

3  k,  Jr 


)• 


In  addition  we  use  the  following  identity,  which  can  be  computed  by  almost  any 
symbolic  system. 


e2xi{^x'+^Xi)dXldx  3 


LiL3  if  pi  =  p3  =  0. 
0  otherwise. 


(17) 


Rather  than  compute  the  projection  of  the  entire  differential  equation,  we  will 
illustrate  the  technique  using  the  following  term  from  (15), 


We  can  ignore  the  Re  u\  term  since  it  is  a  constant. 
Our  goal  is  to  compute  such  that 


CO 


V  Lliz 


2,  E 

1  Jfc,=-oo 


1*3 


(^2)- 


Jfcj  =  -00 


(18) 


The  f[nl  are  functions  of  the  a l")  .  They  are  obtained  bv  taking  inner  products 
(integrals)  of  ( 18)  with  the  orthonormal  basis  functions.  The  first  two  inner  products 
are  Fourier  transforms,  that  is 


ft 


1*3 


e-2xkxxx/ LX 


C-^kyZ./L,^ 


The  final  inner  product  takes  the  form 

/?*,=  /  fkikj 

Jo 

The  final  term  of  U  is  a  polynomial  in  x?  and  is  easily  dealt  with.  The  Fourier 
transform  gives 


-Ilk 

»  *‘.‘j 


Xl\ 


M) 
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(«  (var  X  kl  k3) 

(integral  (*  (sake-saspled-iunctioa 

(laabda  (x)  (-  x  (/  («  x  x)  H)>) 
(dot-product  (eigen  1  kl  k3) 

(conjugate  (eigen  a  kl  k3)))) 

: loser-bound  0 
:upper-bound  X2>) 


Figure  7:  Weyl  code  for  integral 


The  Galerkin  projection  is  just  a  simple  integral,  so 


e*  (Is  - 1)} = “5-tt  l '  (« ■ -  f)  <.,C di 


The  result  many  symbolic  computations  like  this  is  the  system  of  ordinary  dif¬ 
ferential  equations  shown  in  Appendix  A. 


4.3.2  Numerical  Computation 

Having  produced  a  symbolic  system  of  ordinary  differential  equations  like  that 
shown  in  Appendix  A  we  must  still  compute  each  of  the  coefficients.  This  can 
itself  be  a  rather  complex  undertaking  without  the  proper  abstractions.  Consider 
for  instance,  the  a  piece  of  the  sample  term  computed  in  the  last  section: 


<?;- 


(0 


0(n)m 

v**t*3 


dx  2- 


The  product  of  the  two  eigenvectors 


really  means  the  dot  product: 


V'*I*3V‘*1*3  l*l*3VI*l*3  2*1 


(') 


*3  2*  X*3 


4?, 


4”>- . 

>  3»l*3 


Furthermore,  each  of  the  components  of  the  of  the  eigenvectors  are  complex  valued 
functions  that  are  only  known  by  their  numerical  values  at  selected  points. 

To  deal  with  this  problem  we  have  extended  Weyl  [15]  for  this  problem  to  include 
a  new  type  of  a  object,  a  •‘sampled  function.”  A  sampled  function  is  a  function  from 
R  — >  C  (or  R)  that  is  represented  by  its  value  at  certain  points.  When  evaluated  at 
other  points,  its  value  is  automatically  interpolated  (extrapolated)  form  the  values 
at  which  it  is  known.  Like  all  objects  in  Weyl,  arithmetic  with  sampled  functions 
can  be  performed  using  the  usual  Lisp  operators,  including  conjugation. 

The  eigenvectors  <^4"^  are  just  (Weyl)  3-vectors  of  sampled  functions.  Being 
3- vectors  we  can  use  the  dot-product  operator  to  multiply  them.  The  whole  ex¬ 
pression  can  thus  be  computed  as  shown  in  Figure  7.  Notice  that  the  Weyl  code  is 
a  direct  translation  of  the  more  mathematical  form  given  above. 
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The  resulting  system  of  differential  equations  is  somewhat  complex,  as  indicated 
in  Appendix  A.  One  of  these  equations,  when  using  only  a  single  eigenfunction,  has 
the  form 


dj  =  6.1at  +  2.1aia5  4-  l.lajaJ  +  0.4a3a4  4-  O.Sa^aj 

—  (3.0aiaj  4-  3.70203  4-  2.40303  -j-  1.3a4a4  4-  0.60505)01. 


(19) 


where  we  have  only  given  the  coefficients  to  one  decimal  place  for  conciseness.  The 
a,  are  complex  valued  functions,  so  to  numerically  integrate  the  equations,  each 
a,  must  be  converted  to  a  pair  of  real  valued  functions.  For  (19)  this  gives  the 
following  pair  of  equations. 


ii  =  6.1n  -5-  2.1(x2xj  4-  yzyi)  -r  l.l(x3x2  4-  y3yz)  +  0.4(x.sX3  4-  y4y3 j 
t  0.3(x5r4  -f  ysyt)  -  3.0(xf  4-  y\ )x1  4-  3.7(x|  4-  yl)z  1 
~  2.4(x|  4-  y|)x,  4-  1.3(x;  4-  y;)x,  4-  0.6(x|  4-  yf)xIt 
yi  =  6-lyi  -  2.1(x2y1  -  y^x  1)  -  i.i(x3y2  -  jfex2)  -  0.4(x4y3  -  y4x3) 
-  0.3(x5y4  -  y5x4)  -  3.0(x,  4-  y?)yj  4-  3.7(x|  4-  y2)yi 
4-  2.4(x|  4-  y|)yi  4-  1.3(x*  4-  y;)yi  4-  0.6(r|  4-  yf)yi- 


Currently,  we  are  integrating  these  equations  with  the  LSODE  package  {5,  10). 
A  typical  integration  is  shown  in  Figure  8.  There  are  periodic  bursts  of  behavior, 
where  the  equations  become  very  stiff.  We  are  currently  generating  the  Jacobians  of 
these  equations  symbolically  to  speed  the  calculation  during  the  stiff  regimes.  Since 
the  right  hand  sides  of  these  equations  are  polynomials,  symbolic  differentiation 
does  not  cause  the  expressions  to  grew. 

The  bursts  of  activity  in  Figure  8.  when  converted  into  a  velocity  fluctuation, 
correspond  to  the  periodic  formation  of  the  counter-rotating  vortices.  Thus  the 
reduced  system  of  ordinary  differential  equations  has  the  same  qualitative  behavior 
has  the  far  more  complex  Navier-Stokes  equations.  We  are  currently  studying  how¬ 
to  make  this  correlation  more  quantitative  and  how  the  correlation  with  physical 
behavior  depends  upon  the  number  of  empirical  eigenfunctions  used. 

One  should  note  that  for  a  modest  number  of  empirical  eigenfunctions,  the  size  of 
the  system  of  ordinary  differential  equations  becomes  very  large  and  their  formation 
and  manipulation  without  symbolic  techniques  would  be  impractical. 


5  Conclusions 

In  this  paper  we  have  advocated  the  construction  of  special  purpose  simulators  for 
particular  scenes,  rather  than  building  a  general  purpose  simulator.  Towards  this 
end,  we  have  discussed  one  possible  approach  to  the  construction  of  an  environment 
to  enable  the  construction  of  such  simulators.  We  have  particularly  ft  ssed  on  the 
use  of  symbolic  techniques  to  transform  differential  equations  into  executable  code. 
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I 


Figure  8:  Typical  Amplitudes 

We  have  outlined  two  major  areas  in  which  symbolic  computation  can  be  ef¬ 
fectively  used  in  numerical  computations:  (1)  transforming  differential  equations 
into  equations  that  more  accurately  address  the  questions  being  asked  of  the  sys¬ 
tem  under  study,  and  (2)  the  formation  of  the  numerical  integration  code  itself 
from  libraries  of  technique  fragments.  Both  of  these  techniques  suggest  different 
organizations  of  symbolic  computation  systems  than  we  currently  have  available. 
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Abstract 

This  paper  presents  a  framework  for  reasoning  about 
robust  geoemtric  algorithms.  Robustness  is  formally 
defined  and  a  data  structure  called  an  approximate 
polygon  is  introduced  and  used  to  reason  about  poly¬ 
gons  constructed  of  edges  whose  positions  are  uncer¬ 
tain. 

A  robust  algorithm  for  point  location  in  an  approx¬ 
imate  polygon  is  presented.  The  algorithm  uses  only 
the  signature  of  the  point  (not  its  location)  to  de¬ 
termine  whether  the  point  is  inside  or  outside  the 
polygon. 

An  approximate  polygon  could,  by  shifting  its 
edges  back  and  forth  within  their  error  bounds,  in¬ 
duce  a  large  number  of  different  line  arrangements. 
The  cell  Ca  with  signature  a  in  one  such  arrangment 
will  be  different  than  the  cell  C'a  with  signature  a 
in  another  arrangement.  This  paper  proves  that,  re¬ 
gardless  of  their  positions  and  shapes,  the  cells  Ca 
and  C'a  are  always  to  the  same  side  of  the  polygons 
which  induce  their  respective  arrangements. 

Introduction 

Most  geometric  algorithms  assume  that  perfect  “real” 
arithmetic  is  available.  When  these  algorithms  are 
implemented  they  often  fail  because  this  assumption 
is  not  borne  out;  that  is,  these  algorithms  are  not 
robust.  This  failure  occurs  because  either  the  input 
or  the  intermediate  calculations  are  imprecise,  leading 
to  inconsistent  decisions  by  the  algorithm. 

This  paper  presents  a  framework  for  reasoning 
about  robust  geometric  algorithms  which  operate  on 
polygons.  Robustness  is  formally  defined  and  a  data 
structure  called  an  approximate  polygon  is  introduced 
and  used  to  reason  about  polygons  constructed  of 
edges  whose  positions  are  uncertain. 

A  robust  algorithm  for  point  location  in  an  approx¬ 
imate  polygon  is  described.  The  interesting  aspect  of 


this  algorithm  is  that  in  addition  to  the  polygon’s 
position  being  uncertain,  the  point’s  position  in  the 
plane  does  not  have  to  be  known;  only  the  point’s 
signature  is  important  (that  is,  its  left/right  relations 
to  the  edges  of  the  polygon).  The  point  location  al¬ 
gorithm  has  immediate  practical  application  to  solid 
modeling,  particularly  in  the  robust  intersection  of 
polyhedra. 

An  approximate  polygon  could,  by  shifting  its 
edges  back  and  forth  within  their  error  bounds,  in¬ 
duce  a  large  number  of  different  line  arrangements. 
In  each  of  these  arrangements  some  points  with  a 
given  signature  a  may  or  may  not  appear,  and  if 
they  appear,  they  may  be  to  the  interior  or  to  the  ex¬ 
terior  of  the  polygon  which  induces  the  arrangement. 
An  interesting  uniqueness  theorem  is  presented  which 
states  that  in  all  such  line  arrangements,  the  points 
with  signature  a  in  each  arrangement  are  always  to 
the  same  side  of  the  polygon  which  induces  that  ar¬ 
rangement. 

Practical  Applications 

The  point  location  algorithm  has  immediate  practical 
application  to  solid  modeling.  In  particular,  a  solid 
modeler  performing  an  intersection  operation  needs 
to  determine  whether  an  edge  of  one  polyhedron  in¬ 
tersects  a  face  of  another.  This  is  achieved  by  cal¬ 
culating  the  intersection  of  the  edge  with  the  plane 
in  which  the  face  lies,  and  then  asking  whether  this 
point  of  intersection  is  on  the  interior  of  the  polygon 
representing  the  face.  If  the  boundary  of  the  face  and 
the  location  of  the  point  of  intersection  are  known 
precisely  then  this  is  a  trivial  problem. 

However,  polyhedra  usually  have  overconstrained 
faces  and  vertices,  and  the  exact  locations  of  the  ver¬ 
tices,  edges,  and  faces  of  the  polyhedra  can  require 
a  very  large  number  of  bits  to  represent.  Since  the 
input  is  rounded  off  to  a  small  number  of  bits  the  lo¬ 
cations  of  these  features  are  imprecise.  In  addition, 
the  location  of  the  point  of  intersection  can  be  com- 
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pletely  unknown,  particularly  in  ill-conditioned  cases 
in  which  the  intersecting  edge  lies  very  close  to  the 
plane  of  the  face.  Thus  there  is  an  important  practi¬ 
cal  application  for  a  point  location  algorithm  which 
handles  uncertainty  in  the  face  boundary  and  in  the 
point  location. 

An  approximate  polygon  can  represent  a  face 
whose  boundaries  are  not  known  exactly,  and  the 
point  location  algorithm  presented  in  this  paper  can 
determine  whether  a  point  whose  location  is  also  un¬ 
certain  lies  on  the  interior  of  such  a  face.  Since  both 
the  location  of  the  boundary  and  the  location  of  the 
point  are  uncertain,  the  algorithm  must  make  use  of 
some  other  information.  This  other  information  con¬ 
sists  of  the  signature  of  the  point;  that  is,  its  position 
(left  or  right)  with  respect  to  each  edge  of  the  bound¬ 
ary.  It  is  a  reasonable  assumption  that  such  informa¬ 
tion  exists  since  the  signature  is  often  derivable  from 
logical  information  available  in  the  solid  modeler  (for 
example,  see  Karasick’s  modeler  [5]). 

Background 

The  theory  of  approximate  polygons  is  based  upon 
the  “representation  and  model”  approach  of  Hoff¬ 
mann,  Hopcroft,  and  Karasick  [4].  In  this  approach 
the  algorithm  operates  on  a  computer  representation, 
but  presents  output  as  though  it  were  operating  on 
some  mathematical  model  corresponding  to  the  rep¬ 
resentation. 

An  approximate  polygon  is  a  computer  representa¬ 
tion  of  some  real,  mathematical  polygon,  the  model. 
The  model  is  rarely  explicitly  constructed  by  the  algo¬ 
rithm.  An  approximate  polygon  Prep  can  be  thought 
of  as  a  set  of  constraints  on  the  topology  and  position 
of  the  implicit  model  polygon.  Any  real  polygon  P 
satisfying  these  constraints  is  considered  a  model  for 
Prep- 

Under  the  representation  and  model  approach,  the 
definition  of  robustness  is  very  close  to  that  of  Fortune 
(2).  Consider  a  geometric  problem  V  as  a  function 
from  an  input  space  consisting  of  models  to  an  out¬ 
put  space,  V  :2  — ►  O,  and  consider  an  algorithm  A  as 
function  from  a  different  input  space  consisting  of  rep¬ 
resentations  to  the  same  output  space,  A  :  H  — *  O. 
Given  a  representation  xrsp,  the  set  of  its  models  is 
denoted  MODELS(rrep).  This  leads  to  a  definition  of 
robustness: 


An  algorithm  A  for  a  problem  V  is  robust  if 
V  Xrep  €72,  3  x  €  MODELS(xrep) 
such  that  A(Xrep)  =  P(x)- 


Note  that  we  can  pick  an  arbitrary  x  6 
MODELs(xrep).  It  could  be  that  there  are  two  mod¬ 
els  x1  and  x2  such  that  ^(x1)  56  P(x2).  In  this  case 
the  algorithm  could  choose  to  output  either  ^(x1) 
or  V{x2)  and  would  still  be  considered  to  be  robust. 
This  leads  to  a  definition  of  consistency: 


A  problem  V  and  a  representation  7 Z 
are  consistent  if 

V  Xrep  €  72,  V  X1,!2  €  MODELS(Xrep), 

^(x1)  =  7>(x2). 

A  definition  of  correctness  would  be  similar  to  that 
of  robustness,  except  that  the  model  and  representa¬ 
tion  spaces  would  be  identical  and  MOD£LS(xrep)  = 
{Zr-sp}. 

In  evaluating  geometric  algorithms  which  use  the 
representation  and  model  approach,  the  criteria  of 
robustness  and  consistency  should  be  used  in  place  of 
the  usual  criterion  of  correctness. 

Note  that,  unlike  in  Fortune’s  work  [2],  there  is 
no  notion  of  stability  in  the  definition  of  robustness. 
That  is,  there  is  no  notion  of  the  distance  between  the 
representation  xrep  and  the  model  x  which  allows  us 
to  say  “the  implementation  is  stable  if  x  is  near  xrep'’. 
However,  bounds  on  the  models  can  be  achieved  by 
ensuring  that  MODELs(xrep)  is  sufficiently  small. 

Other  approaches 

The  approach  with  approximate  polygons  is  most 
similar  to  that  of  Milenkovic’s  hidden  variable 
method  [6].  His  method  constructs  arrangements  of 
pseudolines  which  are  constrained  to  lie  within  strips 
of  fixed  width.  The  pseudolines  can  be  considered  as  a 
model  and  the  strips  as  a  representation.  Milenkovic’s 
pseudoline  arrangement  algorithm  can  then  be  said 
to  be  provably  robust  in  the  sense  of  the  above  def¬ 
inition.  It  is  interesting  to  note  that,  as  with  many 
algorithms  of  the  “representation  and  model”  variety, 
the  model  is  never  explicitly  constructed. 

There  are  several  other  approaches  to  building  ro¬ 
bust  algorithms  (where  “robust”  is  defined  differ¬ 
ently).  Sugihara  [10,  11,  12]  emphasizes  removing 
redundant  decisions  from  the  algorithm  in  order  to 
maintain  topological  consistency.  Salesin,  Stolfi,  and 
Guibas  [8]  use  what  they  call  epsilon  geometry  to  rea¬ 
son  about  the  amount  of  perturbation  of  the  input  re¬ 
quired  for  certain  epsilon  predicates  to  be  true.  Con¬ 
struction  of  robust  algorithms  is  based  upon  these 
epsilon  predicates.  Dobkin  and  Silver  ( l]  keep  track 
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of  roundoff  error  and,  when  the  error  becomes  too 
large,  increase  precision  and  backtrack  to  some  ear¬ 
lier  point  in  the  computation.  Segal  and  Sequin  [9] 
alter  the  symbolic  data  to  make  it  more  amenable  to 
precise  computation,  and  signal  the  user  when  tol¬ 
erances  on  the  input  become  too  large.  (Milenkovic 
also  alters  the  symbolic  data  in  his  “data  normaliza¬ 
tion”  approach  [6].)  Greene  and  Yao  [3]  discretize  the 
problem  domain,  allowing  the  algorithm  to  perform 
exact  computations. 

In  the  remainder  of  the  paper  approximate  poly¬ 
gons  are  defined,  some  of  their  properties  are  enumer¬ 
ated,  the  point  location  algorithm  is  outlined,  and  the 
uniqueness  theorm  for  point  location  in  an  approxi¬ 
mate  polygon  is  presented. 

Definitions 

A  polygon  P  =  (ej,  ej,  ...e„)  is  an  ordered  list  of  di¬ 
rected  edges,  where  each  edge  e;  lies  on  a  line  £,•  and 
only  intersects  the  edges  e,_i  and  el+i  at  its  end¬ 
points.  Each  edge  is  directed  such  that  the  interior  of 
the  polygon  it  to  its  right. 

An  approximate  polygon  closely  mirrors  the  ap¬ 
pearance  of  a  real  polygon,  as  shown  in  Figure  1. 
The  approximate  polygon  consists  of  an  ordered  list 
of  bands  corresponding  to  the  edges  of  the  model.  The 
position  of  the  bands  in  the  plane  constrains  the  line 
equations  of  the  model. 


Figure  1:  An  Approximate  Polygon 

Just  as  real  polygons  are  based  upon  lines,  approx¬ 
imate  polygons  are  based  on  swaths:  A  swath  5,-  is 
the  region  between  two  lines  If11  and  £jn.  These  lines 
have  the  restriction  that  Vx  The 

restriction  causes  the  lines  to  be  parallel  and  conve¬ 
niently  defines  the  region  between  them  as 

Si  =  {z  |  >  0  A  <r‘(x)  <  0}. 

Just  as  an  edge  is  part  of  a  line,  a  band  is  part 
of  a  swath.  Assuming  for  now  that  an  approximate 


polygon  is  represented  by  an  ordered  list  of  swaths, 
a  band  B{  of  an- approximate  polygon  Prtp  having 
swaths  5,  is  the  shaded  region  in  Figure  2,  defined  as 

Bi  =  Si  n  j  n  £J+i, 

where  F'-  -  /  ^  I  ^  °}  if  '8  COnVeX 

J-  1  {*  |  l)n(z)  >  0}  if  i/j  is  reflex. 


Figure  2:  The  Pieces  of  a  “Band” 

It  will  be  useful  to  define  the  corners  of  a  band 
as  c‘*,  c”,  c*',  and  c*0,  where  h,  t,  and  o  denote 
head,  tail,  in  and  out.  The  definitions  are  shown  in 
the  following  table,  and  depend  on  whether  the  bands 
adjacent  to  the  corner  make  a  convex  or  a  reflex  turn. 
In  Figure  2  the  tail  of  B,  is  convex,  so  the  definitions 
for  c“  and  ct0  are  chosen  from  the  “convex”  column  of 
Table  1.  Since  the  head  of  B,  is  reflex,  the  definitions 
of  c1"  and  c*°  come  from  the  “reflex”  column. 


CONVEX 

REFLEX 

c“ 

n  lft\ 

n  e\h 

ct0 

(cat  n  ^ 

t?'  n  (\h 

c*1' 

3"  n 

C  n  3+i 

cho 

tp*  n  e#\ 

T'  n  3+i 

Table  1:  Defining  the  Band  Corners 

An  approximate  polygon  Prep  is  an  ordered  list  of 
bands1  Bi  which  lie  on  swaths  5;.  Bi  fl  B;-  =  0  iff 
i  and  j  differ  by  more  than  one.  A  real  polygon  P 
is  a  model  for  an  approximate  polygon  Prep  if  the 
following  constraints  are  met: 

1  When  given  a*  input  to  an  algorithm,  the  band*  are  defined 
exactly  with  Boating  point  number*.  Subsequent  computation 
on  the  band*  i*  al*o  done  exactly  (with  extended  precision,  if 
necessary). 
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1.  There  is  a  one-to-one  correspondence  between 
the  lines  U  of  P  and  the  swaths  5,-  of  Prep.  As¬ 
sume  that  corresponds  to  5,-. 

2.  Erich  line  £i(x)  =  0  must  lie  between  the  corners 
of  band  B,-.  That  is,  it  must  satisfy  the  following 
four  constraints  (see  Figure  2): 

£,(c“')<0,  £,(cto)>0, 

ti(chi)  <  0,  £<(cAo)  >  0. 


It  will  be  useful  later  on  to  talk  about  the  span  of 
a  band.  This  is  the  set  of  points  swept  out  by  all 
lines  which  fit  within  the  band.  The  left  and  right 
of  a  band  are  the  set  of  those  points  to  the  left  and 
right  of  the  span.  By  convention,  the  interior  of  the 
approximate  polygon  is  to  the  right  of  the  band.  In 
Figure  3  the  shaded  region  is  span(B,)  and  to  its  left 
and  right  are  LEFT(B;)  and  RIGHT(B,).  For  a  band 
Bi,  define  the  set  of  lines  satisfying  Condition  2  above 
as  LINES(B,). 


SPAN  (Bi) 

right(B,-) 

left(B,-) 


{x  I  3  l  €  LiNEs(Bi),  f(x)  =  0} 
{x  I  V  l  €  UNES(Bi),  t{x)  <  0} 
{x  I  V  l  e  LINES(Bi),  t(x)  >  0} 


Figure  3:  The  span  of  a  Band 


6.  If  x  £  span (B,)  then  in  all  models,  x  lies  to  the 
same  side  of  e,-. 

7.  If  x  £  SPAN(B,+1)  -  B{  then  in  all  models  x  lies 
to  the  same  side  of  e,-. 


Robust  Point  Location  in 
Approximate  Polygons 

The  point  location  problem  would  be  simple  if  the  ex¬ 
act  location  of  the  point  were  given.  However,  in  most 
practical  applications  the  point’s  location  is  known 
only  to  be  within  some  region  of  uncertainty.  In  par¬ 
ticularly  ill-conditioned  situations  this  region  of  un¬ 
certainty  can  be  as  large  as  the  polygon  itself. 

Some  practical  applications  (geometric  modelers, 
for  example)  can,  from  other  information,  logically 
deduce  the  left/right  status  of  the  point  with  re¬ 
spect  to  each  edge  of  the  polygon.  Call  this  l/r 
sequence  the  signature.  If  the  polygon’s  location  is 
known  exactly,  then  in  the  induced  line  arrangement 
a  cell  decomposition  can  easily  determine  whether  all 
points  with  a  given  signature  lie  inside  or  outside  the 
polygon.  It  is  a  different  matter,  however,  when  there 
is  uncertainty  in  the  polygon’s  location.  If  uncer¬ 
tainty  is  modeled  with  an  approximate  polygon  then 
the  following  questions  must  be  answered: 

Question  1  (Robustness)  Given  an 
approximate  polygon  Prtp  and  a  signature 
a  €  (L|R)*,  does  PTtp  have  a  model  P  in 
which  the  induced  line  arrangement  con¬ 
tains  a  cell  with  signature  a,  and  is  the  cell 
INSIDE  or  OUTSIDE  the  model  PI 

Question  2  (Consistency)  Consider 
that  an  approximate  polygon  can  have  two 
models,  P1  and  P3,  which  induce  two  dif¬ 
ferent  line  arrangments.  These  two  arrange¬ 
ments  each  contain  a  cell  with  signature  a 
(call  them  Cl  and  C2).  Then  is  it  possible 
that  C 1  is  inside  Pl  and  C3  is  outside  P3? 


Some  useful  properties  follow  from  the  previous 
definitions  (these  are  stated  without  proof). 

1.  An  approximate  polygon  is  closed  and  simple. 

2.  Edge  e,-  lies  completely  within  B,-. 

3.  Edges  e,-  and  e,-+i  intersect  within  B<  D  B,+ 1- 

4.  All  models  of  an  approximate  polygon  are  sim¬ 
ple. 

5.  SPAN(B,)nsPAN(B,+i)  =  B,  D  B,+1. 


If  the  answer  to  Question  2  were  affirmative  then 
the  signature  a  and  the  approximate  polygon  Prtp 
would  not  be  sufficient  information  to  determine 
point  location,  and  the  problem  would  not  be  con¬ 
sistent.  The  Uniqueness  Theorem  which  is  presented 
later  proves  that  this  is  not  the  case. 

Some  final  definitions 

A  signature  ap(v)  is  a  string  in  (l|r)*.  The  signature 
denotes  the  relation  of  the  point  v  to  each  edge  e, 
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of  the  polygon  P.  The  ith  element  of  ap(v)  is  the 
relation  of  the  point  v  to  edge  e,-  of  the  polygon  P. 
af  (v)  =  R  means  that  v  is  to  the  right  of  edge  e,-  in  P 
and  af  (v)  =  L  means  that  v  is  to  the  left  of  edge  e,- 
in  P.  The  superscripts  will  be  dropped  if  the  polygon 
in  question  is  evident. 

Refer  to  Figure  3  for  the  following  definitions.  A 
half-region  is  similar  to  a  half-space,  except  that  it 
has  a  polygonal  boundary.  The  following  half-regions 
Ri  and  Li  consist  of  those  points  which,  in  at  leati 
one  model  P,  are  either  on  e,-  or  to  the  right  or  LEFT 
of  e,-,  respectively,  in  that  model.  Given  some  a,(v), 
the  half-region  ff,-  is  that  region  in  whose  interior  v 
must  lie  if  it  is  to  have  or,(t> )  as  the  i,l>  component  of 
its  signature.  The  interior  of  the  cell  Ca  consists  of 
those  points  which  haye  signature  a  in  at  least  one 
model. 

Ri  =  SPAN  (P,  )  U  RIGHT(Pi) 

Li  =  SPAN(fli)  U  LEFT(Pi) 

Ri  if  a,-  =  R 
Li  if  a,-  =  L 

a 

Ca  =  f]  Hi 

1=1 

The  next  two  lemmas  will  be  used  to  construct  the 
point  location  algorithm.  The  first  lemma  shows  that 
for  each  point  in  Ca  there  exists  some  model  in  which 
the  point  has  signature  a;  the  second  lemma  shows 
how  to  determine  whether  the  point  is  inside  or  out¬ 
side  that  model. 

Lemma  1  (Model  Existence)  Given  an  approxi¬ 
mate  polygon  Prtp  and  a  signature  a,  construct  Ca  as 
described  above.  Then  for  each  point  v  on  the  interior 
ofCa,  there  exists  some  model  P  €  MODELs(Prep)  in 
which  v  has  signature  a. 

Proof  Since  v  €  Ca,  for  each  i,  v  6  Hi  and  there  is 
some  edge  e,-  in  the  band  Bi  which  has  v  to  the  side 
specified  by  or,-.  These  edges  join  to  form  a  model 
polygon  P  in  which  v  has  signature  a.  □ 

Lemma  2  (Point  Location) 

Given  an  approximate  polygon  Pr<„  a  model  poly¬ 
gon  P  €  MODELS(Prtp),  and  a  point  v  which  has  a 
signature  a  with  respect  to  P,  the  following  are  true: 

J.  If  v  is  strictly  to  the  interior  of  Prtp  (that  is,  it 
does  not  lie  on  any  band  Bi)  then  a,  v  INSIDE  P. 

2.  If  v  it  strictly  to  the  exterior  of  Prtp  then 
V  OUTSIDE  P. 

3.  If  v  e  Bi,  but  o  g  Bin,  then  v  inside  P  iff 

Cti  =  R. 


If  v  €  BiHBi+i  and  the  i/i+l  comer  is  convex, 
then  v  INSIDE  P  iff  a{  =  R  and  af+1  =  R. 

5.  If  v  £  Bi  D  Bi+i  and  the  i/i  +  1  comer  is  reflex, 
then  v  INSIDE  P  iff  a,-  =  R  or  ai+i  =  R. 

Proof  In  Figure  4  the  cases  I  through  5  are  demon¬ 
strated  by  the  points  n  through  x$-  O 


Figure  4:  Cases  for  the  Point  Location  Lemma 

Given  the  Model  Existence  Lemma  and  the  Point 
Location  Lemma,  a  point  location  algorithm  can  be 
developed.  This  algorithm  will  construct  the  region 
Ca,  pick  a  point  from  its  interior,  and  apply  the  rules 
of  the  Point  Location  Lemma  to  determine  whether 
the  point  is  INSIDE  or  OUTSIDE  the  model  in  which  it 
has  signature  a.  The  following  Uniqueness  Theorem 
shows  that  if  one  such  point  is  INSIDE  its  model  poly¬ 
gon  then  all  such  points  are  INSIDE  their  respective 
model  polygons  (similarly  for  outside). 

Theorem  1  (Uniqueness)  Given  an  approximate 
polygon  Prtp  and  a  signature  a,  if  for  some  model 
polygon  in  MODELs(Prep)  there  is  a  point  with  signa¬ 
ture  a  which  is  INSIDE  the  polygon,  then,  for  every 
model  polygon,  all  points  which  have  signature  a  with 
respect  to  that  polygon  are  INSIDE  that  polygon  ( sim¬ 
ilarly  for  OUTSIDE ). 

Proof  in  Appendix  A. 

Point  Location  Algorithm 

The  Model  Existence  Lemma,  Point  Location 
Lemma,  and  Uniqueness  Theorem  combine  to  form 
the  point  location  algorithm  shown  in  Figure  5.  Note 
that  the  algorithm  is  quite  simple  and  never  actually 
constructs  the  model  polygon. 

Lemma  3  (Robustness)  The  point  location  algo¬ 
rithm  is  robust. 


o 


1.  Compute  Ca. 

2.  If  Ca  =  0  then  no  model  of  Prep  induces 
a  cell  with  signature  a. 

3.  Pick  a  point  w  on  the  interior  of  Ca. 

4.  Apply  the  Point  Location  Lemma  to  de¬ 
termine  whether  w  is  inside  or  outside 
of  the  models  in  which  it  has  signature 
a. 


Figure  5:  Point  Location  Algorithm 


Proof  This  follows  directly  from  the  Model  Exis¬ 
tence  Lemma  and  the  Point  Location  Lemma.  □ 

Lemma  4  (Consistency)  The  approximate  point 
location  problem  is  consistent. 

Proof  This  follows  directly  from  the  Uniqueness 
Theorem.  □ 

Lemma  5  (Complexity)  The  point  location  algo¬ 
rithm  has  time  complexity  0(n2). 

Proof  Step  1  of  the  algorithm  can  be  accom¬ 
plished  by  computing  the  arrangement  of  half-regions 
in  0(n2)  time.  This  r  done  by  computing  the  ar¬ 
rangement  of  the  3n  lines  which  define  the  n  half¬ 
regions  Hi,  then  joining  adjacent  cells  which  are  sep¬ 
arated  by  a  line  segment  which  is  not  part  of  the 
boundary  of  some  Hi .  Those  cells  separated  by  a  line 
segment  which  is  part  of  the  boundary  of  some  Hi 
will  have  signatures  which  differ  in  a  single  position, 
so  the  signature  of  each  cell  can  be  found  in  constant 
time. 

The  remaining  steps  of  the  algorithm  take  constant 
time.  Step  3  is  easily  accomplished  given  the  convex 
decomposition  of  Ca  which  is  computed  simultane¬ 
ously  with  Ca  itself.  Thus,  the  overall  running  time 
is  0(n2).  □ 

Summary 

Most  geometric  algorithms  are  not  robust ;  they  fail 
due  to  inexact  input  or  with  inexact  intermediate 
computations.  This  paper  has  introduced  (a)  formal 
definitions  of  robustness  and  consistency,  and  (b)  the 
notion  of  an  approximate  polygon,  along  with  several 
of  its  properties.  With  these,  one  can  formally  de¬ 
velop  robust  and  consistent  algorithms  that  deal  with 
inexact  polygons. 


One  such  algorithm  for  point  location  in  an  approx¬ 
imate  polygon  has  been  presented.  The  algorithm  is 
particularly  suited  for  practical  application  in  a  solid 
modeler  because  it  assumes  uncertainty  in  both  the 
polygon  position  and  the  point  position.  The  point 
location  algorithm  has  been  proved  robust,  and  the 
point  location  problem  has  been  shown  to  be  consis¬ 
tent. 
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Appendix  A 

Theorem  1  (Uniqueness) 

Given  an  approximate  polygon  PTtp  and  a  signature 
a,  if  for  some  model  polygon  in  MODELS(Prep)  there  is 
a  point  with  signature  a  which  is  INSIDE  the  polygon, 
then,  for  every  model  polygon,  all  points  which  have 
signature  a  with  respect  to  that  polygon  are  INSIDE 
that  polygon  (similarly  for  OUTSIDER. 

Proof  (by  contradiction ):  Let  al(z)  be  the  signa¬ 
ture  of  point  x  in  model  P*.  Let  ef  be  edge  e,-  in 
model  Pk.  Then  assume  the  following: 

3 Pl,  P 2  6  MODELs(Prep),  3u,  v  6 

U  INSIDE  P1  A  0  OUTSIDE  P 2  A  O^u)  =  a2(o). 

The  theorem  is  proved  by  showing  that  there  is 
some  edge  e*  which  always  separates  u  and  v,  violat¬ 
ing  this  assumption. 

Lemma  6  One  of  u  or  v  must  lie  within  the  bound¬ 
ary  of  P„p. 

Proof  Assume  that  neither  u  nor  t>  is  within  the 
boundary.  Then  by  the  initial  assumption  and  the 
Point  Location  Lemma  u  and  v  lie  on  opposite  sides  of 
the  boundary.  Say  u  is  inside  and  v  is  outside.  Then 
the  segment  uv  must  traverse  both  of  the  parallel 
sides  of  some  band  B„  as  shown  in  Figure  6.  From 
the  definition  of  the  corners  of  B,  and  the  definition 
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of  span  (Bi),  u  and  v  lie  to  different  sides  of  span  (B,). 
Then  by  Property  6  u  lies  to  one  side  of  all  models  of 
Prep  and  v  lies  to  the  other  side.  Thus,  in  the  models 
P1  and  P3,  a}(u)  ^  aj(v)  and  the  initial  assumption 
is  contradicted.  O 


v 


Figure  6:  uv  crosses  some  B,- 


From  Lemma  6  assume  without  loss  of  generality 
that  «  lies  in  the  boundary  of  Prtp. 

Lemma  7  The  point  u  lies  in  some  band  Bt-  such 
that  a}(u)  =  R  and  ctf(v)  =  R. 

Proof  u  is  INSIDE  Pl  by  the  initial  assumption 
and  is  in  some  band  by  Lemma  6.  If  u  6  B*  and 
u  £  Bt±  j,  then  by  the  Point  Location  Lemma  u  is  to 
the  right  of  e£.  Choose  i  =  k.  If  u  €  Bt  flflt+i  then 
by  the  Point  Location  Lemma  u  is  to  the  right  of  at 
least  one  of  ej  or  e£+1.  Choose  i  to  be  k  or  k  +  1  to 
satisfy  the  lemma.  Then  by  the  initial  assumption, 
a^u)  =  R  means  that  a?(v)  =  R  also.  0 

Lemma  8  On  the  segment  W  there  is  some  point 
z  £  v  which  is  inside  P3.  Furthermore,  e3  does  not 
intersect  uv  between  z  and  v. 

Proof 

Case  1:  u  6  B,-,  but  u  £  B,± j.  If  a?(u)  =  R, 
then  u  INSIDE  P3  (by  the  Point  Location  Lemma),  so 
choose  x  —  u  and  we  are  done.  Otherwise  consider 
the  case  in  which  a3(u)  ±  R. 

Refer  to  Figure  7.  The  point  u  lies  in  B,-  and  not  in 
Bj_  j,  so  by  Property  7,  in  all  models  Pk,  u  is  to  the 
same  side  of  ef_l.  In  particular,  a}_ t(u)  =  aj_j(u). 
By  the  initial  assumption,  q}_i(u)  =  of_i(v),  so 
Q?-i(u)  =  By  a  similar  argument  for  e,+i, 


Figure  7:  x  exists  if  u  €  B, 

Ql+i(u)  —  Q?+i(°)>  Thus  u  and  v  lie  between  two 
lines  touching  the  endpoints  of  e?. 

Since  af(u)  £  r  and  since  by  Lemma  7  af(v)  =  R, 
uv  must  cross  the  line  defined  by  e?.  But  u  and  v  lie 
between  the  two  lines  touching  the  endpoints  of  e3, 
so  W  must  cross  the  edge  ef .  Since  by  Property  4  all 
models  P  are  simple,  there  is  a  small  neighborhood 
to  the  right  of  e }  which  contains  only  points  interior 
to  P.  The  segment  uF  passes  through  this  neighbor¬ 
hood,  so  there  is  some  point  x  jb  v  on  uv  which  is 
INSIDE  P3,  and  Ttmcf  =  0. 


Figure  8:  z  exists  if  u  €  B,-  f)  B,+i 

Case  i:  u  €  B,nBi+1.  Refer  to  Figure  8.  By  an 
argument  similar  to  Case  1,  u  and  r  must  lie  to  the 
same  side  of  ef_j  in  all  models  Pk,  and  must  lie  to 
the  same  side  of  ef+3  in  all  models  Pk.  If  u  inside  P 3 
then  choose  x  =  u  and  we  are  done.  Otherwise,  if  u 


OUTSIDE  P2  then  the  edges  e2  and  e2+1  must  separate 
u  and  u.  Thus  the  segment  uv  must  cross  either  ej  or 
e2+1,  and  in  doing  so  must  pass  through  a  neighbor¬ 
hood  of  interior  points  to  the  right  of  the  edge  that 
it  Ciosses.  Therefore  there  is  some  point  x  ^  v  on  ut; 
which  is  inside  P2,  and  xv  n  e2  =  0.  □ 

Lemma  9  /n  P2  there  is  some  edge  ej  which  crosses 
uv  such  that  aj(u )  =  R  and  ctj(v)  =  L.  Furthermore, 
ej  can  be  chosen  such  that  no  other  ej  crosses  uv 
between  ej  and  v. 

Proof  From  lemma  8,  uv  contains  a  point  x  which 
is  inside  P2,  and  by  the  initial  assumption,  v  out¬ 
side  P2.  From  the  Jordan  curve  theorem  we  know 
that  xv  must  cross  the  boundary  of  P2  on  some  edge 
ej  with  x  to  the  inside  (right)  of  e2  and  v  to  the  out¬ 
side  (left)  of  ej.  From  the  ordering  of  points  along 
uv  (u  <  x  <  v),  if  aj(x)  =  R  then  cej(u)  =  R  also.  If 
there  are  many  candidates  for  ej,  choose  that  which 
is  closest  to  t;  to  satisfy  the  second  part  of  the  lemma. 
□ 

Lemma  10  Bj  n  £,•  =  0 

Proof  If  |i  —  j|  >  1  then  by  the  definition  of  an 
approximate  polygon  Bi  n  Bj  =  0.  So  we  only  have 
to  consider  |i  —  j\  <  1.  But  by  Lemma  9  ej  intersects 
uv  between  x  and  v,  and  by  Lemma  8  ej  does  not 
intersect  uv  between  x  and  v.  So  ej  ^  ej  and  hence 
*  5*  h 

Assume  that  j  —  i—  l.  Then,  by  the  initial  assump¬ 
tion,  aj(v )  =  L  means  that  orj(u)  =  l.  By  Lemma  9, 
a2(u)  =  R.  Since  u  is  to  different  sides  of  e,  in  P1 
and  P2,  u  must  lie  in  span  (Bj).  Since  u  also  lies  in 
Bi,  by  Property  5  u  lies  in  the  corner  B n  Bj. 

Suppose  corner  i/j  is  convex.  Since  u  INSIDE  P1,  by 
the  Point  Location  Lemmatt-(«)  =  Rand  aj(u)  =  R. 
But,  by  the  initial  assumption,  aj(u)  =  R  means  that 
aj(v)  =  R,  contradicting  lemma  9.  So  corner  i/j  is 
not  convex. 

Suppose  corner  i/j  is  reflex.  Refer  to  Figure  9.  By 
Lemma9,  aj(u)  =  Rand  crj(v)  =  l,  and  by  Lemma 7, 
aj(v )  =  R.  But  if  ut?  i*  to  intersect  ej  then  it  must 
also  intersect  ej  closer  to  v,  aa  shown  in  the  Figure. 
This  contradicts  Lemma  9,  which  states  the  ej  is  the 
closest  intersection  to  v.  So  corner  i/j  is  not  reflex. 

Thus  the  assumption  is  false  that  j  jt  i  —  l.  By  a 
similar  argument  j  £  i  +  I.  Therefore  Bj  n  Bi  =  0. 
□ 

Lemma  11  u  €  SPAN  (Bj)  -  Bj 

Proof:  By  lemma  9,  aj(u)  =  R  and  aj(v)  =  L.  By 
the  initial  assumption, aj{v)  =  L  means  that  oj(u)  = 


Figure  9:  Utf  must  intersect  ej  closer  to  v 


L.  Since  ory(u)  is  different  in  models  Pl  and  P2.  u 
must  lie  in  span(£;).  From  lemma  10,  u  cannot  lie 
in  Bj  since  it  already  lies  in  □ 

Lemma  12  Define  k  such  that  Bj  n  Bt  is  closest  to 
u.  Then  u  and  v  are  on  opposite  sides  of  ej. 

Proof  Refer  to  Figure  10.  Note  that  k  =  j  ±  1, 
otherwise  Bj  and  Bj,  wouldn’t  intersect  at  all.  In  any 
model  P  the  point  u  and  the  edge  e;  are  on  opposite 
sides  of  the  line  defined  by  et  (u  £  span  (Bj)-  B,  by 
lemma  11,  and  since  Bk  is  closest  to  u,  it  separates 
u  from  the  rest  of  Bj,  which  contains  ej).  Thus  uv 
must  intersect  ej  at  some  point  z  to  the  side  of  e\ 
which  is  opposite  to  u.  From  the  ordering  of  points 
along  uv  ((u  <  z  <  u),  v  must  also  be  opposite  to  u. 
□ 


Lemma  13  orj(u)  =  ctj(u). 

Proof  By  Lemma  12,  u  £  SPAN  (Bj)  —  Bj.  By 
Property  5,  SPAN(£j)nsPAN(f?*)  =  BtflBi.so  with 
a  bit  of  algebra  we  can  conclude  that  u  $  SPAN(Pt). 
Then  by  Property  6,  aj(u)  =  a2(u).  □ 


By  lemma  12  there  is  some  edge  e2  in  model  P2 
such  that  o|(u)  ^  a 2(t>),  and  by  lemma  13,  oj(u)  = 
qJ(u).  So  a£(u)  £  a\{ v).  But  this  contradicts  the 
initial  assumption,  so  the  theorem  is  proved  by  con¬ 
tradiction. 

O 
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Figure  10:  tt  and  v  are  on  oppoeite  sides  of  ej 
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Abstract 

This  paper  presents  a  polynomial  time  algorithm  for  determining  whether  a  given  univariate 
rational  function  over  an  arbitrary  field  is  the  composition  of  two  rational  functions  over  that 
field,  and  finds  them  if  so. 


1  Introduction 

The  problem  of  determining  if  a  function  can  be  written  as  the  composition  of  two  “smaller’’  functions 
/( i)  =  g(h(x))  has  been  of  interest  for  a  long  time.  Until  now,  work  has  focused  on  the  univariate, 
polynomial  version  of  this  problem:  When  can  the  polynomial  f(x)  be  written  as  g(h(x)),  where 
both  g(x)  and  h(x)  are  polynomials?  The  original  work  in  the  symbolic  computation  community 
was  presented  in  1976  [2],  but  the  algorithms,  which  in  the  worst  case  required  exponential  time, 
were  not  published  until  1985  [3],  This  was  soon  followed  by  the  work  of  Kozen  and  Landau  [1 1]  who 
provided  a  polynomial  time  algorithm  for  decomposition  of  polynomials  over  fields  of  characteristic 
zero,  which  did  not  require  factorization  o,  polynomials.  Some  additional  improvements  and  analysis 
of  the  positive  characteristic  case  where  then  presented  by  von  zur  Gathen  [23,  21,  22].  A  number 
of  other  papers  have  since  been  published  on  different  extensions  and  variations  of  this  problem  [1, 
7,  5,  4]. 

All  of  these  results  deal  with  polynomial  decomposition.  The  generalization  to  rational  functions, 
which  has  significantly  wider  applicability,  appears  to  be  a  far  harder  problem.  Notice  that  in  the 
polynomial  case,  the  degree  of  g  and  h  must  divide  the  degree  of  /.  This  limits  the  number  of 
different  polynomials  that  must  be  considered  and  even  allows  one  to  solve  the  problem  by  looking 
for  solutions  of  non-linear  algebraic  equations  (admittedly  in  exponential  time).  When  /,  g  and  h 
are  rational  functions,  there  is  no  immediately  obvious  bound  on  the  degrees  of  the  numerators  of 
g  and  h ,  since  the  numerator  and  denominator  of  g(h(x))  could  have  a  common  factor.  In  fact,  no 
such  common  factor  can  arise,  as  we  prove  below. 

Furthermore,  we  demonstrate  that  in  the  rational  function  case,  g  and  h  can  be  determined  from 
/  in  polynomial  time.  This  algorithm  is  valid  even  if  the  characteristic  of  the  field  is  positive, 
which  for  the  polynomial  case  is  not  a  completely  resolved  problem.  One  other  difference  between 
our  approach  and  other  approaches,  is  that  in  this  paper  we  obtain  a  decomposition  over  the  fie:d 
of  definition  of  f(x).  Thus  we  may  fail  to  find  rational  function  decompositions  that  exist  over 
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algebraic  extensions.  Such  issues  do  not  arise  for  the  corresponding  problem  of  polynomials  over  a 
field  of  characteristic  zero,  but  do  for  polynomials  over  fields  of  finite  characteristic. 

Section  2  provides  some  general  background  material.  In  Section  3  we  present  the  new  algo¬ 
rithms  for  rational  function  decomposition.  Finally,  we  comment  on  previous  work  in  and  give  some 
conclusions  in  Section  4. 

2  Preliminaries 

Let  f(x)  be  a  rational  function  in  x  with  coefficients  in  a  field  k.  We  extend  the  notion  of  degree  of  a 
polynomial  by  defining  the  degree  of  f(x),  denoted  by  deg  /,  to  be  the  maximum  of  the  polynomial 
degrees  of  the  (relatively  prime)  numerator  and  denominator  of  /.  The  degree  of  the  field  k(x)  over 
k(f(x))  is  the  degree  of  /,  if  /  is  a  polynomial.  This  remains  true  even  if  /  is  a  rational  function, 
as  shown  by  the  following  proposition. 

Proposition  1  Let  k(x)  be  an  extension  of  the  field  k(f(x))  where  f(x)  is  a  rational  function  of 
degree  n.  Then  [k(x):k(f(x))]  =  n. 

Proof:  Denote  the  numerator  of  f(x)  by  p(x)  and  the  denominator  by  q(x).  We  can  instead  consider 
the  isomorphic  fields  k(y)  —  k(f(x))  and 

*(y)W/(p(z)  “  </?(*))  -  *(*)• 

P(x,  y)  =  p(x )  —  yq(x)  is  primitive  as  a  polynomial  in  y  since  p(x)  and  q(x)  are  relatively  prime. 
Since  it  is  linear  in  y  it  is  irreducible.  Therefore,  the  degree  of  x  over  the  field  k(y)  is 

degr  P(x, y)  =  max(degp, deg?)  =  deg/. 

□ 

Let  f(x)  =  g{h(x))  be  a  rational  function  decomposition  over  a  field  k.  The  following  proposition 
provides  bounds  on  the  degrees  of  g(x)  and  h(x)  in  terms  of  the  degree  of  f{x).  In  principle,  this 
result  gives  an  algorithm  for  rational  function  decomposition,  albeit  an  exponential  time  algorithm. 

Proposition  2  Assume  f(x),  g(x)  and  h(x)  are  elements  of  k(x)  such  that  f(x)  =  g(h(x)).  Then 

deg  /  =  (deg  g)  •  (deg/i) 

Proof: 

Consider  the  fields  shown  in  Figure  1.  The  degrees  of  the  extensions  are  [k(x):k(h(x))]  =  deg/t, 
[Jt(x):  Jt(/(x))]  =  deg  /  and  [ifc(y)  :k (y(y))]  =  degy.  k(h(x))  is  an  algebraic  extension  of  k(f(x)) 
inside  k(x).  Thus, 

deg  /  =  (*(r)  :*(/(*))] 

*[*(*):*(*(«)))  *[4(M*))  :*(/(*))] 

=  [k(x):k(h(x))}.[k(y):k(g(y))) 

=  (degh)  •  (degy), 

using  Proposition  1.  D 

A  function  that  is  the  ratio  of  to  linear  polynomials  is  called  a  fractional  linear  function,  viz. 

A(x)  =  (or  +  b)/(cx  +  d). 


k(x) 

k(Kx))-*— —  k(y) 

*(/(*))- —  *(?(y)) 

Figure  1:  Fields  involved  in  decomposition 


Fractional  linear  functions  have  degree  1.  If  two  fields  k(fi(x))  and  k(fi(x))  are  isomorphic  then 
there  exist  rational  functions  such  that 


/i(*)  =  fli(/2(*)) 

/*(*)  =  /M/i(*))  =  fl2(/?i(/2(x))) 


By  Proposition  2  (deg  Ri)  ■  (deg  R2)  =  1  and  R\  and  Rj  must  be  fractional  linear  functions. 

We  say  that  two  rational  functions  are  linearly  equivalent  if  there  exists  fractional  linear  functions 
Ai  and  A2  such  that 

f(x)  =  \1(g(X2(x))). 

Two  decompositions  (polynomial  or  rational  function) 


/  =  ffi  0  92  o  •  •  * 0  9m 
=  hi  o  h2  o  •  •  •  o  h„ 


are  said  to  be  equivalent  if  m  =  n  and  <7,-  is  linearly  equivalent  to  hi. 

The  link  between  field  structure  and  function  decomposition  comes  from  Luroth 's  theorem ,  which 
was  proven  by  Luroth  [15]  for  k  =  C  and  by  Steinitz  in  general  [18], 

Proposition  3  (Luroth)  //  k  £  K  C  k(x)  then  K  =  k(g(x))  where  g(x)  is  a  rational  function  in 
x  over  k. 


An  elementary  proof  of  Luroth’s  theorem  may  be  found  in  van  der  Waerden  [20].  An  effective 
proof  appears  in  Weber  [24]  §124,  and  in  English  in  Schinzel  [17]. 

The  key  insight  in  studying  functional  decomposition  is  the  following  corollary  of  Luroth’s  the¬ 
orem. 


Proposition  4  Let  k  be  an  arbitrary  field  and  f{x)  a  rational  function  over  k.  There  is  a  one  to  one 
correspondence  between  the  lattice  of  subfields  between  k(x)  and  k{f(x))  and  the  rational  function 
decompositions  of  f(x)  up  to  equivalence. 

Proof:  If  f(x)  has  a  nontrivial  decomposition  f(x)  =  g(h(x))f  then  k(h{x))  will  be  an  intermediate 
field  between  k(x)  and  k(f(x)).  Conversely,  if  K  is  field  intermediate  between  k(x)  and  k(f(x)) 
then  it  must  be  of  the  form  k(h(x)),  where  h(x)  is  a  rational  function  in  x.  k(h(x))  is  canonically 
isomorphic  to  k(y)  as  shown  in  Figure  1,  where  <Ph(y)  *-*  h(x).  k(f(x))  is  intermediate  between 


k(y)  =  k(h(x))  and  k,  so  by  Luroth's  theorem,  there  is  a  rational  function  g(y)  such  that  k{f(x))  = 
k(g(y)).  Thus  f(x)  is  linearly  equivalent  to  g(h(x)).  D 

The  following  two  propositions  follow  from  Proposition  2  and  are  quite  useful. 

Propositions  Let  k  be  an  arbitrary  field  and  g  i  and  <72  relatively  prime  elements  of  F[x].  Then 
for  all  polynomials  h(x)  £  k[x],  gi(h(x))  and.  g2(h(x))  are  relatively  prime. 

Proof:  Without  loss  of  generality  assume  that  degyi  >  deg  g.  Define  g(x)  to  be  the  ratio  of  y^x) 
and  ga(x).  Since  yt  and  gi  are  relatively  prime  and  deg <71  >  deg  <72,  degy(x)  =  degyi.  Let 


f(x)  =g(h(x))  = 


9i(h(x))  fi(x) 


92(h(x))  h(x)' 

where  f\  and  f%  are  relatively  prime.  Thus 

deg  fi(x)  <  degy,(/i(x))  =  (deg?,)  •  (deg  /i), 

where  equality  holds  if  and  only  if  yi(h(x))  and  S2(/»(x))  are  relatively  prime.  Furthermore,  deg/i  > 
deg/2  so  deg  /  =  deg/i-  By  Proposition  2 

deg  /(x)  =  (deg  g)  ■  (deg  h)  =  (deg  y  1)  ■  (deg  h) 

so  deg/i(x)  =  (degyi)  •  (deg/i)  and  gi(h(x))  and  y2(/»(x))  are  relatively  prime.  D 

The  argument  of  previous  proposition  applies  equally  when  h(x)  is  a  rational  function.  In  this 
case,  it  is  best  to  view  g^  and  <72  as  bivariate  homogeneous  functions  of  the  same  degree,  which  gives 
the  following  result. 

Proposition  6  Lei  g  1  and  <72  be  relatively  prime,  homogeneous  polynomials  in  two  variables.  If  hi 
and  are  also  relatively  prime  polynomials,  then  0i(/»i,/i2)  and  02(hi,/»2)  are  also  relatively  prime. 

Notice  that  the  requirement  that  g  1  and  <72  be  homogeneous  is  necessary  as  the  following  example 
shows: 

5i(x,y)  =  x+  1  ’j 
9a{x,y)  =  y-2 

MO  =  t 
MO  =  + 1 J 

As  a  consequence  of  Proposition  6,  rational  function  decomposition  can  be  viewed  as  a  coupled 
polynomial  decomposition  problem,  viz. 

/i(*,y)  =  0i(Mz.y).M*.y))> 
h{i,y)  =  y2(Mz,y),Mz,y)), 

where  /,,  g,  and  h,  are  homogeneous  polynomials  and  the  pairs  {/i,/2},  {01,02}  and  {/»i,/»2}  have 
the  same  degree. 


0  i(MM  —  t  +  1 

02(/»j,/»2)  =  t2-  1 


3  Rational  Function  Decomposition 

The  bounds  of  Proposition  2  provide  significant  insight  into  rational  function  decomposition.  In 
particular,  if  the  degree  of  /(x)  is  prime,  then  it  has  no  non-trivial  decomposition.  A  simple,  expo¬ 
nential  time  algorithm  for  determining  a  decomposition  can  be  constructed  by  using  undetermined 


coefficients.  Assume  that  deg  /  =  rs  and  we  are  looking  for  a  decomposition  f(x)  =  g(h(x)),  where 
deg  <7  =  r  and  deg/i  =  s.  We  can  write  g  and  h  in  terms  undetermined  coefficients,  e.g. 

g(x)  =  9n =  goxr  +gixr~l  +  ...  +  9r 

9d(x)  gr+lXr  +gr+2Xr~l  +  •••  +  <72r+r 

There  are  2r  +  2  undetermined  coefficients  in  g(x)  and  2s  +  2  in  h(x).  By  Proposition  6,  we  can 
treat  the  numerator  and  denominator  of  f(x)  independently.  Equating  the  coefficients  of  x'  in  the 
following  equations  gives  a  system  of  2 rs  +  2  algebraic  equations  in  the  <?,■  and  ht. 

fox"  +•••  +  /„ 

=  9ohn(x)r  +---  +  grhd(x)r 

fr,  +  lXr>  H - +  fjrt+l 

=  gr+ihn(x)r  +  •  •  •  +  92r+iMx)r 

Any  decomposition  of  f(x)  is  a  solution  to  this  system  of  equations.  Conversely,  any  solution  to  this 
system  for  which  deg$r  =  r  and  deg/i  =  s  gives  a  decomposition  of  f(x).  However,  this  approach  is 
not  very  efficient.  Nonetheless,  it  does  demonstrate  the  existence  of  an  algorithm. 

The  efficient  techniques  that  have  been  developed  all  tend  to  be  divided  into  two  phases,  com¬ 
puting  h(x)  and  then  given  h(x)  computing  g(x).  (The  hard  part  is  finding  h(x).)  We  discuss  the 
phases  out  of  order  for  simplicity.  Determining  g  from  /  and  h  is  discussed  in  Section  3.1,  while  the 
determination  of  h  is  discussed  in  Section  3.2. 

3.1  Determining  g  from  /  and  h 

The  most  direct  way  to  obtain  g(x)  such  that  /(x)  =  g(h(x)),  when  /  and  h  are  known  is  to  explicitly 
solve  the  linear  equations  for  the  coefficients  of  j(x)  that  arise  from  (1).  This  approach  is  discussed 
in  detail  by  Dickerson  [5,  4]  as  “computing  the  left  composition  factor.”  In  this  section  we  present 
a  simple  analytic  technique  that  relies  on  reversion  of  power  series  and  is  valid  when  the  coefficient 
field  has  characteristic  0. 

Let  Xj  be  a  fractional  linear  function  such  that  /  =  A/  o  /  has  a  zero  at  0.  Define  h  and  A/, 
similarly.  If  /  =  g  o  h  then 

/(x)  =  (AJ1  ogoXh)oh{x), 

and  g(x)  =  (A^ 1  ojo  A>,)(r).  So  without  loss  of  generality  we  can  assume  /( 0)  =  /i(0)  =  0. 
h(x)  has  a  power  series  expansion  of  the  form 

h(x)  =  htxl  +  ht+ix1*1  +  •  •  • 

Using  standard  techniques  [10]  we  can  obtain  a  power  series  in  t  for  x  in  t  =  h(x) 

x  =  h-l(t)  =  h\t1/l  +  h,2tvt  +  ---. 

Replacing  x  by  this  power  series  in  the  power  series  for  /(x)  we  get  a  power  series  in  t.  If  there 
are  any  fractional  powers  then  there  does  not  exist  a  “left  composition  factor.”  Compute  the  first 
2r  terms  of  the  power  series  expansion  of  /(/»-1(x))  at  0.  The  (r,  r)  Pade  approximate  [16]  to  this 
power  series  is  the  only  possible  candidate  for  g(x).  This  power  series  technique  may  be  easier  to 
program  than  Dickerson’s  technique,  and  using  fast  power  series  techniques  [12]  it  might  have  better 
asymptotic  complexity. 


*(*) 


E[a]/(f(a)-t)  =  F 


k(h(X)) - Emm-t) 


*(/(*)) - m  =  e 


Figure  2:  Field  Structure 


3.2  Determination  of  h 


For  rational  function  decomposition,  we  determine  h(x)  by  explicitly  determining  a  subfield  of  k(x) 
and  then  use  a  constructive  version  of  Luroth’s  theorem  to  compute  a  generator  for  the  subfield. 
The  tower  of  fields  we  will  be  working  with  is  shown  in  Figure  2.  Note  that  the  fields  on  the  same 
horizontal  line  in  Figure  2  are  isomorphic.  By  Proposition  3  every  subfield  of  F  is  of  the  form 
k(h(x))  and  there  exists  a  rational  function  g  such  that  g(h(x))  =  f(x),  since  k(f(x))  lies  between 
k(y)  =  k(h(x))  and  k.  Thus  every  non-trivial  subfield  of  F  yields  a  non-trivial  decomposition  of 
/(*)• 

To  illustrate  our  procedure  consider  the  following  example: 


/(*)  = 


( +  M 

W-2J  °  \x*  +  2j 


2xA  +  6x2  +  5  _  fn(x) 
z4  -f-  6x2  +  7  fd{x) ' 


where  /„  and  fj  are  relatively  prime.  We  want  to  find  an  intermediate  field  between  Jfc(r)  and 
k(f(x)).  Our  first  step  is  to  convert  these  fields  to  a  more  conventional  form.  If  E  =  k(t)  =  k(f(x)) 
and  £7[or]  =  k(x )  then  a  satisfies  the  minimal  polynomial 


f{t,  Z )  =  fn(Z )  -  tfd(Z)  =  (t  +  2 )Z*  +  (61  +  6 )Z3  +  It  +  5. 


This  polynomial’s  factorization  over  £(a]  is 

f(t,  Z)  =  (Z-  a)(Z  +  a)((t  +  2 )Z2  +  (t+  2)q2  +  6(t  +  1)).  (2) 

Over  a  proper  subfield  of  E[a],  f(t ,  Z)  will  not  factor  so  much.  In  particular,  over  a  subfield  it 
cannot  have  a  linear  factor.  Given  (2),  the  only  poesible  factors  of  /(/,  Z)  over  the  subfield  £T[/3]  are 
Z  -  a1  and  (( t  +  2 )Z‘i  +  (t  +  2)o2  +  6 (t  +  1)).  Thus  E[0]  must  contain  the  coefficients  of  these  two 
polynomials.  If  E{j3 j  is  the  smallest  subfield  of  E[a]  for  which  /(f,  Z)  has  such  a  factorization,  then 
it  must  be  generated  by  the  coefficients  of  these  two  polynomials.  In  this  case  we  can  assume  that 
(3  =  a2,  whose  minimal  polynomial  is 

h(t,Z)  =  (t+2)Zi +  (6t  +  %)Z +  7t  +  S.  (3) 

To  convert  E\fl\  back  to  the  form  k(f(x))  we  observe  that  the  elements  of  E[0\  are  rational 
functions  in  x  over  k  by  Luroth’s  theorem.  When  t  is  replaced  by  f(x),  (3)  must  have  linear  factors, 
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i(/(x),z)  =  (z-*i)(z-gii), 

which  leads  to  the  intermediate  fields  k(x2)  and  k(( 3z2  +  4)/(2z2  +  3)).  These  two  fields  are 
isomorphic  by  the  fractional  linear  map  z  t—  (3z  +  4)/(2z  +  3).  Using  the  k(x2)  as  the  intermediate 
field,  we  have  h(x)  =  x2,  and  thus  the  irreducible  decomposition: 

2z *  +  6z2  +  5  _  _2z2  +  6z  +  5  2 

z*  +  6z2  +  7  z2  +  6z  +  7 

The  original  decomposition  is  equivalent  to  this  one  since 


z2  +  1  /  2z2  +  6z  +  5\  /  — 2z  +  1  \ 
z2  —  2  \  x2  +  6x  +  7  )  \  i— 1  j 

This  basic  approach  is  applicable  to  the  general  problem  except  for  deciding  which  factors  of 
f(t,  Z )  should  be  recombined  to  generate  a  factorization  over  a  subfield  of  ITfor].  We  could  try  all 
possible  combinations  of  factors  of  f(t,  Z)  until  we  find  one  that  yields  a  proper  subfield  of  £”[a]. 
However,  in  the  worst  case  this  would  require  an  exponential  number  of  trials.  Instead,  we  use  a 
version  of  Landau  and  Miller’s  algorithm  BLOCKS  in  [14]  to  find  a  non-trivial  block,  which  will 
generate  a  proper  subfield  of  £[<*].  As  pointed  out  by  Kozen  and  Landau  [11],  this  algorithm  only 
requires  that  the  extension  FfaJ/E  be  separable.  Kozen  and  Landau  may  need  to  examine  as  many 
as  0(nlogn)  non-trivial  blocks  to  find  a  decomposition.  However,  in  our  case,  any  non-trivial  block 
will  give  a  rational  function  decomposition.  These  techniques  allow  us  to  decide  which  factors  of 
f{t,  Z)  should  be  recombined  in  polynomial  time. 

Furthermore,  observe  that  Trager’s  polynomial  time  reduction  of  factorization  over  algebraic 
extensions  [19],  which  was  used  by  Landau  to  show  that  factoring  over  algebraic  number  fields  is 
polynomial  time  [13]  is  applicable  here  also,  so  the  factorization  of  f(t,  Z)  over  the  function  field 
£[«]  can  be  done  in  polynomial  time. 

The  coefficients  of  such  a  factorization  generate  the  intermediate  field  E[0].  Since  we  are  seek¬ 
ing  any  intermediate  field,  a  single  coefficient  that  is  not  in  E  suffices.  The  minimal  polynomial 
of  for  that  coefficient  can  be  determined  using  resultants  and  square  free  decompositions  to  give 
E{j3]/(pff{t,P)).  h(x)  is  then  deduced  from  a  linear  factor  of  pp(/(z),Z),  which  need  only  be  fac¬ 
tored  over  k.  (Factoring  bivariate  polynomials  is  polynomial  time  by  Kaltofen  [9].) 

It  is  worth  commenting  on  the  practicality  of  this  algorithm.  Its  dominant  cost  is  the  factorization 
of  f(t,  Z)  over  £(<)[<*],  which  is  about  as  costly  as  factoring  a  polynomial  of  degree  (deg/)2.  Given 
the  practical  difficulties  of  factoring  polynomials  of  degree  greater  than  about  100,  it  seems  that  it 
will  be  very  difficult  to  determine  the  decomposition  of  /(z)  if  the  degree  of  /(z)  is  greater  than 
about  10. 

3.3  Characteristic  p  case 

Determining  any  decomposition,  as  opposed  to  determining  a  decomposition  with  a  particular  degree 
pattern  over  a  field  of  characteristic  p  is  only  slightly  more  difficult  than  the  characteristic  0  case, 
using  the  technique  of  Section  3.2.  Assume  that  charF  =  p  and  /(z)  is  a  rational  function  over  k. 
The  decomposition  of  f(x)  may  no  longer  be  unique,  but  Proposition  4  shows  that  there  is  still  a  one 
to  one  correspondence  between  the  inequivalent  decompositions  of  f(x)  and  the  fields  intermediate 
between  k(x)  and  i(/(z)). 


k(x) 


So  =  *(/(*)) 


Figure  3:  Field  Structure  for  f(x)  =  xp3+p3  -  xp3+l  -  xp3+p  +  xp+1 


Referring  to  Figure  2,  let  f(t,  Z)  be  the  (irreducible)  minimal  polynomial  of  a  over  E.  If  f(t,  Z) 
is  separable,  then  £[a]  is  separable  over  E  and  a  subfield  can  be  computed  using  the  techniques  of 
the  previous  section.  If  f(t,  Z)  is  inseparable  then  it  can  be  written  as 

/(t,Z)  =  /(*,Zp'‘), 

for  some  positive  value  of  fi.  Furthermore,  /  is  separable  over  E.  Clearly,  the  field  £[ap“j  lies 
between  £[a]  and  £  and  thus  a  linear  factor  of  /(/(x),  Z)  will  give  a  decomposition  factor  of  /(x). 
Since  £[apM]  is  separable  over  £,  the  techniques  of  the  previous  section  can  be  used  to  find  additional 
right  decomposition  factors.  Left  decompositions  factors  can  be  found  from  the  fields  £[ap],  which 
lie  between  £(a]  and  £[ap*]  for  1  <  i  <  p. 

It  is  worth  noting  that  even  the  pathological  example  suggested  by  Dorey  and  Whaples  [6] 

/(x)  =  xp+1  o  (xp  +  x)  o  (xp  -  x), 

=  (xp3  —  Xp3_p+l  —  xp  +  x)  o  xp+1, 

=  xp5+p3  -  xp3+1  -  xp3+p  +  xp+1, 
is  can  be  handled  straightforwardly,  since  the  derived  polyncmial 

f(t,  z)  =  zp5+p3  -  zp3+p  -  zpl+1  +  Zp+1  -  t 

is  separable.  The  fields  associated  with  the  two  decompositions  of  /(x)  are  shown  in  Figure  3. 

In  the  case  of  polynomial  decomposition,  notice  that  /(f,  Z)  is  inseparable  if  and  only  if  /(x)  is 
a  rational  function  of  xp.  Thus  the  distinction  made  by  von  zur  Gathen  [21,  22]  between  “tame” 
and  “wild"  might  more  appropriately  be  made  on  whether  or  not  /(x)  is  a  rational  function  in  xp. 

Note  that  this  approach  only  finds  some  decomposition  of  /(x).  It  cannot  find  a  prescribed  one. 
In  particular,  if  one  is  looking  for  a  decomposition  /(x)  =  g(h(x))  where  p|  deg  g  then  the  extension 
k(x)/k(f(x))  may  be  inseparable  and  we  would  thus  Lave  no  algorithm  for  finding  intermediate 
fields.  This  problem  is  raised  in  [22]. 


4  Conclusions 


The  technique  is  used  to  find  the  h(x)  in  Section  3.2  is  reminiscent  of  the  technique  proposed  by 
Kozen  and  Landau  [11]  for  decomposition  over  arbitrary  fields.  However,  they  studied  intermediate 
fields  between  k(a)/(f(a))  and  k.  While  there  is  an  intermediate  field  between  k(a)  and  ife  whenever 
f(x)  is  decomposable,  the  existence  of  an  intermediate  field  does  not  guarantee  a  decomposition.  By 
using  intermediate  fields  between  fields  jfc(t)[a]/(/(a)-t)  and  k(t),  we  avoid  much  of  the  complexity 
of  their  approach  since  any  such  intermediate  field  does  lead  to  decomposition  of  f(x). 

It  is  tempting  to  conjecture  that  Propositions  5  and  6  can  be  generalized  to  more  variables,  but 
the  straightforward  generalization  is  not  true,  as  pointed  out  in  Section  2.  It  would  be  interesting 
to  know  in  what  way  it  can  be  generalized. 

This  work  has  benefited  from  discussions  with  Barry  Trager  and  Dexter  Kozen.  Susan  Landau’s 
comments  on  an  earlier  version  of  this  paper  where  quite  helpful.  The  diagrams  in  this  paper  were 
typeset  using  Paul  Taylor’s  commutative  diagram  macros  for  DTgX. 
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Abstract 

Modeling  physical  objects  with  low-degree  algebraic  surfaces  shows  promise  for  applications 
where  manipulating  and  reasoning  about  physical  objects  are  important.  In  this  paper,  we 
present  an  algorithm  for  free-form  surface  constructions  using  implicitly  defined  cubic  surface 
patches.  The  input  data  for  the  algorithm  is  an  arbitrary  polyhedron  with  a  normal  prescribed 
at  each  vertex  of  the  polyhedron.  Using  a  Clough-Tocher  like  splitting  scheme,  the  algorithm 
constructs  a  smooth  piecewise  cubic  surface  interpolating  the  vertices  of  the  polyhedron  and  the 
prescribed  normal  at  each  vertex.  The  free-form  surface  construction  in  the  algorithm  is  local 
and  quadratically  precise.  In  addition,  the  shape  of  the  free-form  surfaces  can  be  manipulated 
through  a  set  of  intuitive  shape  parameters  without  knowing  the  details  of  the  algorithm.  The 
implementation  results  are  reported. 

Keywords:  Geometric  modeling,  object  representation,  free-form  surface,  Bernstein-Bezier 
representation,  implicit  patch,  design. 

1  Introduction 

While  developing  a  geometric  modeling  system  for  representing,  manipulating  and  reasoning  about 
physical  objects,  we  derived  and  implemented  an  algorithm  for  constructing  geometric  models  for 
smooth  objects  of  arbitrary  shapes  and  topologies.  Such  geometric  models  are  important  for  solid 
modeling,  computer-aided  design,  visualization,  computer  graphics,  and  robotics. 

The  geometric  models  of  arbitrary  smooth  objects  are  represented  by  closed  free-form  surfaces. 
The  algorithm  we  drive  generates  a  free-form  surface  from  the  input  data  of  an  arbitrary  polyhedron 
with  a  normal  prescribed  at  each  vertex  of  the  polyhedron.  Using  a  Clough-Tocher  like  splitting 
scheme,  the  algorithm  constructs  a  smooth  piecewise  cubic  surface  interpolating  the  vertices  of  the 
polyhedron  and  the  prescribed  normals. 

The  algorithm  we  derive  has  the  following  features.  First,  the  algorithm  is  local,  so  modifying  a 
piece  of  input  data  affects  only  nearby  surface  patches.  Second,  the  algorithm  has  quadratic  preci¬ 
sion,  which  means  that  if  the  input  data  is  taken  from  a  quadric  surface,  the  algorithm  reproduces 
the  quadric  surface.  Finally,  the  shape  of  the  free-form  surfaces  produced  by  the  algorithm  can  be 
controlled  through  a  set  of  intuitive  shape  parameters  without  knowing  the  details  of  the  algorithm. 

An  important  motivation  of  our  w’ork  is  to  construct  geometric  models  that  facilitate  manipulat¬ 
ing  and  reasoning  about  physical  objects  (Hopciuft  and  Ivrafft  1986;  Hoffmann  1989).  Traditionally, 
the  building  blocks  for  free-form  surface  constructions  are  parametric  patches.  As  far  as  design  and 
display  are  concerned,  parametric  patches  are  very  successful.  But  when  it  comes  to  manipulating 
and  reasoning  about  physical  objects,  parametric  patches  run  into  serious  problems.  Parametric 
patches  are  not  closed  under  some  elementary  operations  in  geometric  modeling,  such  as  sweeping 
and  Minkowski  sum  (Bajaj  and  Kim  1987).  The  intersection  of  two  parametric  patches  is  extremely 
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difficult  to  represent  and  evaluate  (Hoffmann  1980)  because  the  algebraic  degree  of  the  intersection 
is  prohibitively  high.  As  an  example,  we  notice  that  in  general  the  intersection  of  two  commonly 
used  bicubic  patches  is  a  space  curve  of  degree  324  (Hopcroft  and  Krafft  1986). 

These  problems  can  be  avoided  by  building  free-form  surfaces  from  low-degree  implicit  patches. 
Implicit  patches  are  closed  under  all  common  operations  required  by  a  geometric  modeling  system 
(Bajaj  1989),  and  the  intersection  of  two  degree  n  implicit  patches  has  degree  n2,  which  is  small  if 
n  is.  Low- degree  implicit  patches  also  allow  the  use  of  algebraic  techniques  as  opposed  to  numerical 
techniques  in  reasoning  about  physical  objects  (Hopcroft  and  Krafft  1986).  These  features  make 
implicit  patches  a  superior  choice  for  applications  where  manipulating  and  reasoning  about  physical 
objects  are  important.  In  addition,  from  a  practical  point  of  view,  implicit  patches  are  compact  to 
store  and  relatively  easy  to  ray  trace. 

An  inviting  class  of  implicit  patches  for  free-form  surface  constructions  is  the  class  of  quadric 
patches.  When  the  input  data  is  a  polyhedron  without  normals  prescribed  at  its  vertices,  a  free¬ 
form  surface  can  be  constructed  using  quadric  patches.  However,  quadric  patches  have  fundamental 
limitations  that  make  it  impossible  to  allow  prescribing  normals  in  the  input  data.  Roughly  speak¬ 
ing,  when  a  free-form  surface  is  constructed  by  replacing  the  facets  of  the  input  polyhedron  with 
quadric  patches,  they  introduce  a  correlation  between  the  normals  at  the  adjacent  vertices  of  the 
input  polyhedron.  We  have  investigated  the  role  of  quadric  patches  as  primitives  for  free-form 
surface  constructions,  and  we  hope  to  report  the  results  elsewhere. 

Being  able  to  prescribe  the  normals  in  the  input  data  is  important.  Prescribing  normals  is  a 
measure  to  control  the  patches  in  the  free-form  surfaces  so  that  only  a  few  patches  are  needed  for 
representing  a  smooth  object  that  would  otherwise  requires  thousands  of  polygons  to  approximate. 
One  way  to  overcome  the  limitations  of  quadric  patches  is  to  split  the  edges  of  the  input  polyhedron, 
as  was  done  by  Dahmen  (Dahmen  1989).  However,  from  a  theoretical  point  of  view,  Dahmen's 
method  cannot  handle  arbitrary  input  polyhedron  because  his  method  requires  the  existence  of 
“transversal  systems”,  which  no  one  know:,  how  to  construct  in  general;  from  a  practical  point 
of  view,  splitting  the  edge  of  the  input  poljliedron  causes  oscillations  in  the  free-form  surfaces, 
making  it  impossible  to  produce  free-form  surfaces  of  pleasing  shapes.  In  this  paper,  we  show  that 
the  limitations  of  quadric  patches  can  be  overcome  by  cubic  patches. 

1.1  Previous  work 

There  is  a  rich  literature  on  surface  constructions  using  parametric  patches,  and  a  recent  survey  can 
be  found  in  (Mann  et  al.  1990).  Modeling  complex  objects  with  implicit  patches  was  introduced 
in  recent  years  and  is  becoming  an  increasingly  prominent  area  of  research.  General  techniques 
for  implicit  modeling  are  developed  by  researchers  vorldwide  (Xishimura  et  al.  1985;  Bloomenthal 
and  Wyvill  1990;  Dahmen  1989),  In  particular,  many  authors  have  demonstrated  the  power  of 
implicit  patches  in  deriving  blending  surfaces  (Blinn  1982;  Middleditch  and  Sears  1985;  Hoffmann 
and  Hopcroft  1987;  Rockwood  and  Owen  1987)  and  in  surface  fitting  and  approximation  (Bajaj 
and  Ihm  1989;  Patrikalakis  and  Kriezis  1989). 

Sederberg  proposed  using  Bernstein-Bezier  representation  of  implicit  patches  in  free-form  sur¬ 
face  constructions  (Sederberg  1985).  Subsequently,  various  techniques  are  developed  for  construct 
ing  free-form  surfaces  using  implicit  patches  ( Patrikalakis  and  Kriezis  1989;  Bajaj  and  Ihm  19S9; 
Moore  and  Warren  1990;  Sederberg  1990).  In  particular,  Patrikalakis  et  al..  Sederberg.  and  Bajaj 
et  al.  demonstrated  the  complications  and  pitfalls  of  modeling  with  implicit  patches  (Patrikalakis 
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and  Kriezis  1989;  Bajaj  and  Him  1989;  Sederberg  1990). 

Dahmen  (Dahmen  1989)  gave  an  algorithm  for  constructing  free-form  surfaces  from  quadric 
patches.  But  the  algorithm  cannot  handle  arbitrary  polyhedra,  and  the  splitting  scheme  in  the 
algorithm  prevents  it  producing  pleasing  shapes.  There  are  also  algorithms  for  constructing  free¬ 
form  surfaces  with  implicit  patches  of  degree  six  and  degree  five  (Moore  and  Warren  1990;  Bajaj 
1990). 

2  Conceptual  overview 

The  algorithm  described  in  this  paper  builds  free-form  surfaces  from  the  input  data  of  a  polyhedron 
with  a  normal  vector  prescribed  at  each  vertex  of  the  polyhedron.  The  input  data  is  denoted  by 
(V,  A'”),  where  V  is  an  arbitrary  polyhedron  with  vertex  set  {xj,  ■  •  -  and  X  is  a  set  of  normals 
{ni,  •  •  • ,  n;t}  with  n,  being  the  normal  vector  prescribed  at  x,.  The  facets  of  V  are  assumed  to  be 
triangular. 

The  basic  idea  of  the  algorithm  is  very  simple.  The  free-form  surface  to  be  built  must  be 
in  the  neighborhood  of  the  input  polyhedron  V,  so  we  construct  a  neighborhood  S  of  V  using 
tetrahedra  and  creat  a  cubic  polynomial  for  each  tetrahedron  used.  By  ensuring  C 1  conditions 
between  adjacent  tetrahedra,  we  obtain  a  global  C1  function  that  is  a  cubic  polynomial  in  each 
tetrahedron.  The  zero  contour  of  this  global  C1  function  within  the  neighborhood  E  is  the  free-form 
surface  to  be  generated. 

The  following  three  aspects  are  crucial  to  the  success  of  the  algorithm. 

1.  The  construction  of  a  neighborhood  E  of  the  input  polyhedron  using  tetrahedra.  The  neigh¬ 
borhood  must  locally  contain  the  tangent  plane  determined  by  the  prescribed  normal  at  each 
vertex  of  the  input  polyhedron  V ,  and  the  neighborhood  must  have  the  same  topology  as  the 
polyhedron  V. 

2.  A  scheme  for  defining  a  globally  Cl  function  which  is  a  cubic  polynomial  over  each  tetrahedron 
within  the  neighborhood  S.  The  scheme  must  leave  free  control  points  in  the  definition  of 
each  cubic  polynomial  so  that  the  zero  contour  of  the  cubic  polynomial  can  be  controlled  by 
these  free  control  points. 

3.  A  mechanism  to  control  the  cubic  polynomial  defined  for  each  tetrahedron  so  that  the  zero 
contour  of  the  cubic  polynomial  inside  the  tetrahedron  is  a  single-sheeted  cubic  patch  without 
holes,  extraneous  sheets,  self  intersections,  or  other  topological  anomalies. 

These  three  aspects  will  be  stressed  throughout  the  development  of  the  algorithm. 

3  Algorithm  details 

Now  we  address  the  three  aspects  of  the  algorithm  in  detail.  In  this  paper,  we  use  [xi---x„]  to 
denote  the  convex  hull  of  point  set  {xj,  •  •  -  .x,,}. 

3.1  The  construction  of  the  neighborhood  E 

The  basic  spatial  elements  used  to  build  the  neighborhood  E  of  the  polyhedron  V  are  tetrahedra. 
Tetrahedra  are  chosen  for  two  reasons,  one.  tctiahedra  are  simple  and  flexible  three  dimensional 
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Figure  1:  Filling  gaps  between  two  double  tetrahedra 


spar  3  units;  two,  tetrahedra  facilitate  the  use  of  Bernstein-Bezier  representation,  which  is  the  base 
for  this  work. 

The  neighborhood  E  is  constructed  as  follows.  For  each  facet  F  =  [X1X2X3]  of  the  polyhedron 
V,  two  _  ints  X4  and  y4  off  each  side  of  the  facet  are  chosen,  and  they  determine  two  tetrahe¬ 
dra,  [X1X2X3X4]  and  [x1x2X3y4].  These  two  tetrahedra  form  a  double  tetrahedron  denoted  by 
([X1X2X3X4],  [x1x2x3y4]).  Consider  an  adjacen*-  facet  F'  =  [X1X2X3]  and  its  double  tetrahedron 
([xix2x3X4],  [x,1X2X3y^]).  Between  the  double  tetrahedra  of  facets  F  and  F',  there  are  two  gaps. 
One  gap  is  between  the  tetrahedra  [XJX2X3X4]  and  [X1X2X3X4]:  the  other  is  between  [xiX2X3y4] 
and  [xjX2X3y^].  The  first  gap  is  filled  with  a  pair  of  tetrahedra  [x"x2X3X4]  and  [XJX2X3X4],  and 
the  second  gap  is  filled  with  another  pair  of  tetrahedra.  [y"x2X3y4]  and  [y"x2X3y^].  Here  x"  and 
Yi  are  points  on  the  line  segments  [X4X4]  and  [y^y^ j  respectively.  All  these  are  shown  in  Figure  1. 

As  an  auxiliary  geometric  structure  for  the  free  form  surface  construction,  the  neighborhood  E 
must  satisfy  the  following  condition.  At  each  vertex  x,  of  the  polyhedron  V,  the  neighborhood  E 
should  locally  contain  the  tangent  plane  defined  by  n,.  In  other  words,  there  is  a  disk  D  around 
the  vertex  x,-  in  the  tangent  plane  at  x,-  such  that 

DCS. 

3.2  A  scheme  for  enforcing  C1  conditions  over  S 

Having  built  a  neighborhood  E  of  the  polyhedron  V.  we  construct  a  C1  function  /  over  the  neigh¬ 
borhood  E  so  that 

/(x,-)  =  0,  V/(x,}  =  n,-,  i  =  1, •••,£.  (1) 

The  zero  contour  of  /  within  E  is  the  free-form  surface  to  be  generated. 

The  construction  of  /  can  be  outlined  as  follows.  First,  we  split  the  the  tetrahedra  that  have 
facets  of  V  as  faces:  the  neighborhood  S  is  kept  the  same  except  some  of  its  tetrahedra  are  split. 
Then,  the  function  /  is  defined  by  constructing  a  cubic  polynomial  for  each  tetrahedron  within  the 
neighborhood  E. 

To  show  the  splitting  scheme,  we  take  a  facet  [X|X>X3]  and  its  double  tetrahedron  ([xiX2X3X4j. 
[xiX2X3y.j])  as  an  example.  Let  w  be  a  point  in  the  facet  [X]X>X3].  We  split  the  double  tetrahedron 
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□  control  points  to  be  decided 

O  control  points  that  are  free  or  are  decided  by  free  control  point. 

■  control  points  from  the  input  data 
®  control  points  from  C1  conditions 

Figure  2:  The  C1  conditions  between  two  adjacent  double  tetrahedra 

into  six  tetiahedra:  [x,xJwx4]  and  [x,x^wy4]  for  1  <  i  <  j  <  3.  For  symmetry  and  robustness 
reasons,  w  is  often  chosen  to  be  the  centroid  of  triangle  [X1X2X3],  while  y^  and  X4  are  chosen  to  be 
on  the  line  that  passes  through  w  and  is  perpendicular  to  [X1X2X3]. 

The  construction  of  cubic  pov  normals  over  the  tetrahedra  within  S  takes  two  steps.  Consider  a 
facet  F'  —  [xiX2X3]  adjacent  to  facet  F  -  [xix^]  and  the  double  tetrahedron  of  F',  ([X1X2X3X4], 
[xjxaxay^]).  The  facet  F'  and  its  double  tetrahedron  are  split  at  the  centroid  w'  of  F'  in  the 
same  way  that  ([X1X2X3X4],  [xiX2X;j.y.i])  is  split  at  w.  For  the  facets  F  and  F\  the  first  step  takes 
place  over  tetrahedra  Vt  sr  [X2X3X4W],  V2  =  [x2X3X4w'| ,  Wj  =  (x2X3x/,,X4],  W2  =  [x2X3x"x!!,], 
K  =  [x2x3y4wj,  V{  =  fx2xcy^w'j,  W{  =  [x2x3yi'y-{]»  and  W'2  =  [x2x3y"y4]  as  in  Figure  2.  We 
construct  the  cubic  polynomials  over  tetrahedra  lTt,  If  2,  W[ ,  and  ftj.  At  the  same  time,  the  cubic 
polynomial^  over  tetrahedra  \2,  V{,  and  1 '2'  are  partially  determined  through  C1  conditions, 
The  same  process  is  carried  out  between  every  pait  of  adjacent  facets  of  V,  so  at  the  end  of  the  first 
step,  the  cubic  polynomials  over  the  tetrahedi  a  [x,x_,X4w]  and  [x,Xjy4w]  are  partially  constructed 
for  all  i  <  j  <  3.  l»ieu,  the  second  step  completes  the  construction  of  these  cubic  polynomials 


according  to  C1  conditions. 

Now  we  describe  the  first  step  in  detail.  Throughout  this  description,  we  assume  i  =  1,2 
whenever  i  appears.  By  doing  so  we  are  taking  advantage  of  the  symmetry  of  the  problem  in 
consideration. 

Let  the  cubic  polynomials  /,  over  V,,  /,'  over  Vj',  g,  over  W,,  and  g[  over  \V't  be  expressed  in 
Bernstein-Bezier  forms  as  follows. 


fi(x)  =  J2  WliTi), 

|A|=3 

(2) 

9i(x)  =  E  Wl (W). 

|A|=3 

(3) 

//(*)=  E  c\Bl(r[), 

|A|=3 

(4) 

g\{x)  =  E  ABlipW 
|A|=3 

(5) 

where  r,,  t(,  p[  and  p,  are  the  barycentric  coordinates  on  Vj,  Vf,  and  W,  respectively.  We  call 
the  ax’s,  b\  s,  ca’s,  and  d\s  the  control  points  of  the  cubic  polynomials  /,,  /,',  g„  and  g[  respectively. 
Our  task  is  to  determine  these  control  points. 

For  notaticnal  convenience,  if  two  tetrahedra  sharing  a  common  face,  we  equal  the  control 
points  of  the  associated  cubic  polynomials  on  the  common  face  to  ensure  C°  continuity.  Hence 
such  control  points  will  be  defined  only  once. 

All  the  control  points  over  tetrahedra  V,,  F/,  IF,,  and  \V[  that  can  be  determined  from  the 
input  data  are  as  follows.  The  fact  that  the  zero  contours  of  /,,  /',  g, ,  and  g[  pass  through  X2  and 
X3  implies 

a0300  =  a0030  = 
c0300  =  c0030  =  Oi 
^0300  =  ^0030  =  0, 


and 


‘0300 


‘0030 


=  0. 


More  control  points  are  determined  by  the  normals  at  the  vertices  X2  and  X3.  For  example, 


aL+e‘  =  "  xj)>  J  =  2»3- 

Similar  expressions  are  used  to  determine  the  control  points  a^e^+e'1  ^or  i  =  2i3  an(^  k  —  1, 4,  ce J+c4 
for  j  =  2,3,  b'2eJ+e ,  for  j  =  2,3,  and  d'2ej+ei  for  j  =  2,3. 

Before  determining  the  rest  of  the  control  points  according  to  C1  conditions,  we  have  to  choose 
some  control  points  to  be  free  control  points  whose  values  will  be  left  unspecified  at  this  point.  This 
is  because  creating  a  piecewise  cubic  C1  function  /  over  the  neighborhood  S  is  only  an  intermediate 
step  in  the  free-form  surface  construction.  Having  defined  a  cubic  polynomial  whose  zero  contour 
passing  through  some  vertices  of  a  tetrahedion  does  not  guarantee  the  existence  of  a  taut  cubic 


G 


Figure  3:  A  handle  on  a  cubic  patch 

patch  inside  the  tetrahedron.  There  may  not  be  a  cubic  patch  inside  the  tetrahedron  at  all,  or  even 
there  is,  the  cubic  patch  can  have  self  intersections,  holes,  and  extra  sheets.  Figure  3  is  a  more 
dramatical  example:  a  handle  appears  on  an  otherwise  nice  cubic  patch.  If  this  cubic  patch  is  in  a 
free- form  surface  constructed  from  the  input  polyhedron  V )  the  topology  of  the  free-form  surface 
is  bound  to  be  different  from  that  of  V. 

We  choose  the  following  free  control  points  a'2et+eJ  {j  =  1,2, 3, 4),  c^4+eJ  (j  =  1,2, 3, 4),  b'200V 
and  c^ooi  f°r  /«,  f',t  <7«>  and  g[  respectively.  The  intuition  of  these  control  points  are  as  follows. 
Control  points  al2e<+eJ  {j  —  1, 2, 3)  are  equivalent  to  the  function  values  and  gradients  at  x*  and  xj, , 
and  the  control  point  62001  enables  us  to  have  complete  control  of  the  function  values  of  g,  along  the 
line  segment  [x4x^].  The  same  statement  can  be  made  about  control  points  c'2et+e)  ( j  =  1,2, 3, 4) 
and  ^2001-  3.3,  we  will  explain  how  these  fiee  control  points  affects  the  associated  cubic  patches. 

Now  we  determine  the  rest  of  the  control  poiats  to  ensure  C1  conditions.  Consider  the  C1 
conditions  across  faces  [X2X3X.1]  and  [x2X3xf,].  Suppose 


x"  =  /3}x  1  +  3\x2  +  A3X3  +  fi\  X'i 


i 


and 


x"  =  P\ fx'i  4-  /?2x2  4-  ^3X3  4-  /?4x!,. 
Then,  the  C1  conditions  are  the  following. 


6*1002  -  / 

'la*1002  +  r^2a0\02  +  /^3a0012  +  ^4a0003) 

^'llOl  =  / 

^a'llOl  +  02aO2Ol  +  $3a0111  +  /^4a0102) 

(6) 

^'loii  =  / 

^ia'l0n  +  $2abui  +  /^3a0021  +  $4a0012) 

(7) 

^*1110  =  / 

^Ia'lll0  +  p2a02l0  +  /^3a0120  +  04aOlll- 

(8) 

can  be  viewed  as  the  definitions  for  the  control  points  6'n01, 

6joiii  and 

and 


6*mo,  leaving  a\011  and  a’1101  to  be  determined.  Equation  (8)  will  be  treated  later. 

Moving  on  to  the  C1  conditions  across  [X2X3X/|/],  we  see  that  if 

X'/  =  /ilX.i  +  H2  x'4, 

then  the  C1  conditions  are  the  following. 

63000  =  162001  +  ^2^2001  > 

^1020  =  Mi^lioi  +  fobiioi, 
b‘l200  —  +  ^2^  1011) 

and 

&mo  =  Miami  +/i2flom 

Again,  the  first  three  equations  can  be  viewed  as  definitions  for  control  points  630001  wio20> 


(9) 

(10) 

(11) 


(12) 
6’i020>  an^ 

6)200!  an(l  the  last  equation  will  be  treated  later.  Notice  that  6\0ii  and  6'1101  in  the  above  equations 
are  defined  earlier  by  the  equations  (6)  and  (7). 

Finally,  we  consider  the  C1  conditions  across  faces  [x2X3y4],  [x2X3y"],  and  [x2X3y)].  All  the 
conirol  points  of  o',  and  some  of  the  control  points  of  /,'  can  be  fixed  in  the  same  way  as  the  control 
points  of  fi  and  <7,-.  In  doing  so,  we  also  have  two  equations  left  untreated. 


^illO  -  7l«ill0  +  72a0210  +  73a0120  +  74c0111 


and 


^'lllO  _  C01 1 1  +  ^Coili) 

where  the  coefficients  t/’s  and  7’s  come  from  the  following  relations: 


(13) 

(14) 


y"  =  7?Xi  4-  7|x2  +  7^X3  +  7,^4 


and 

y"  =  /‘iy-i  +  M2Y4- 

Now  we  collectively  treat  the  equations  (8),  (12),  (13),  and  (14)  as  promised.  These  equations 
can  be  rewritten  as 


(15) 
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and 


Mlfloill  +  ^2001 1 1  — /^l^'llio  + /^2a0210  + /^3a0120  + /^4a0111‘  (16) 

Here  CqU1  ca,n  be  determined  from  a" s  through  the  Cl  conditions  across  [X1X2X3]  and  [X1X2X3], 

cOttt  =  fttatllO  +  Q2a0210  +  C*3a0120  +  a4a0111>  (1^) 


where  the  a’s  come  from 

y4  =  a|xi  +  Q2X2  +  03X3  +  <*4X4 

and 

y4  =  afxi  +  a|x2  +  a^x3  +  0^X4. 

The  equations  (15),  (16),  and  (17)  form  a  system  of  six  linear  equations  with  six  unknowns, 
aom>  amo>  an<*  coxir  When  the  points  X4,  X4,  y4,  and  y^  are  in  general  position,  the  system 
always  has  a  solution.  It  may  happen  that  theie  is  a  family  of  solutions,  in  which  case  we  choose  a 
solution  as  follows.  Using  the  degree-elevation  pioperty  of  Bernstein-Bezier  representation,  we  can 
compute  default  values  for  a'1110  and  CqU1  from  the  prescribed  normals.  A  solution  to  the  system 
can  be  selected  from  the  family  of  solutions  according  to  these  default  values. 

This  finishes  the  first  step  in  the  construction  of  the  cubic  polynomials  over  the  tetrahedra 
within  the  neighborhood  E.  Now  g,  and  g[  are  completely  constructed,  while  /,  and  /'  are  partially 
constructed. 

The  second  step  completes  the  construction  of  /,  and  f[.  For  this  purpose,  we  shift  our  focus  to 
the  double  tetrahedron  ([X1X2X3X4],  [xiX2X3y4]),  which  has  been  split  into  tetrahedra  [X1X2X4W], 
[X1X3X4W],  [X3X2X4W],  [x3x2y4w],  [x3x2y4w],  and  [x3x2y4w]. 

Consider  the  problem  of  completing  the  construction  of  the  partially  constructed  cubic  poly¬ 
nomials  for  tetrahedra  U\  =  [X2X3X4W],  U2  =  [X1X3X4W],  and  U3  =  [xjX2X4w].  Denote  the  cubic 
polynomials  for  Z7,-  by 

hi{vi)  =  Y,  «!\ 

1*1-3 

where  ul  is  the  barycentric  coordinate  of  L\ .  It  is  easy  to  recognize  that  polynomial  hi  is  the  same 
as  fi  in  the  first  step.  More  generally,  the  partially  constructed  functions  h,  is  the  result  of  carrying 
out  the  first  step  for  [xiX2X3]  and  the  facet  sharing  edge  [xmx„]  (m  i,  n  ^  i,l  <  m,n  <  3)  with 
[x  1X2X3].  We  denote  f\,  X'j,  and  T\  by  h\ ,  U\.  and  u\  in  the  second  step  to  reflect  the  new  symmetry. 

The  tasl.  of  ensuring  C1  conditions  between  V,  \s  is  greatly  simplified  by  taking  advantage  of 
the  fact  that  w  6  [xiX2X3j.  The  control  points  over  L\  can  be  divided  into  four  groups.  The  i-tli 
group,  called  i-th  layer,  is  the  set  of  su°h  ^iat  ^4  =  z-  Because  w  6  [xiX2X3],  the  C1 

conditions  between  U,  and  U3  only  involve  contiol  points  from  the  same  layer.  So  we  can  satisfy 
the  C1  conditions  by  examining  each  layer  as  if  we  were  working  on  bivariate  polynomials. 

For  the  0-th  layer,  the  control  points  A  Aj0  are  defined  previously  for  all  Aj  <  1.  Determining 
the  rest  of  the  control  points  in  this  layer  is  exactly  the  famous  Clough-Tocher  interpolation  in 
finite  element  analysis.  Figure  4  illustrates  a  standard  solution  (Farin  1986). 

For  the  I-tli  layer,  the  control  points  «A  A  A  ,  are  defined  earlier  for  all  Ai  =  0.  Since  this  layer 
can  be  viewed  as  a  bivariate  quadratic  function,  the  known  control  points  uniquely  determine  the 
rest  of  the  control  points  within  the  layer  through  the  Cl  conditions  (Farin  1986). 

The  control  points  in  the  2-th  and  3-tli  layers  ate  trivially  determined  by  the  the  function  value 
and  gradient  at  x.j. 


9 


o  centroid  of  surrounding  boxes 
□  constructed  from  previous  step 
v  centroid  of  surrounding  circles 

Figure  *  The  Clough- Tocher  bivariate  splines 

To  complete  the  second  step,  we  carry  out  the  same  argument  for  tetrahedra  [xiX2y4w], 
[x2x3y4w],  and  [xiX3y4w].  As  for  the  C1  conditions  across  [X1X2X3],  notice  that  these  condi¬ 
tions  only  involve  control  points  from  the  0-th  layer  and  the  1-th  layers.  From  equation  (17)  and 
the  way  the  control  points  in  the  1-th  layer  aie  determined,  it  is  easy  to  see  that  the  C1  conditions 
across  [X1X2X3]  are  indeed  satisfied. 

Therefore,  we  have  constructed  the  global  C1  function  /  satisfying  (1).  If  the  free  control  points 
are  chosen  so  that  a  "nice”  cubic  patch  is  obtained  inside  each  tetrahedron  within  the  neighborhood 
S,  then  the  zero  contour  of  /  inside  the  neighborhood  E  is  the  free- form  surface  to  be  constructed. 

3.3  Obtaining  and  controlling  the  cubic  patches 

As  we  mentioned  earlier,  creating  a  Cx  function  /  over  the  neighborhood  S  according  to  (1)  is  only 
an  intermediate  step.  In  general,  such  a  function  rarely  yields  the  free-form  surface  we  expect.  The 
problem  is  that  some  of  the  control  points  of  a  cubic  polynomial  strongly  affect  the  zero  contour 
of  the  cubic  polynomial  inside  the  associated  tetrahedron.  If  we  let  these  control  points  be  decided 
by  the  C1  conditions,  then  the  zero  contour  of  the  cubic  polynomial  inside  the  tetrahedron  exhibits 
various  behaviors  undesirable  for  free-form  surface  constructions. 

The  following  situations  may  occur  for  the  zero  contour  of  a  cubic  polynomial  inside  a  tetrahe¬ 
dron. 

1.  There  is  no  zero  contour  in  th  ’n  trior  of  the  tetrahedron  even  though  the  zero  contour  is 
known  to  pass  through  several  vertices  of  the  tetrahedron. 

2.  There  are  self-intersection  points,  or  singular  points  on  a  cubic  patch. 

3.  There  are  holes  on  a  cubic  patch  caused  by  the  zero  contour  of  the  cubic  polynomial  leaving 
and  coming  back  to  the  tetrahedron.  See  the  left  figure  in  Figure  5. 

4.  There  are  multiple  sheets  of  the  zero  contour  inside  the  tetrahedron.  See  the  left  figure  in 
Figure  6. 
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Figure  5:  Avoiding  holes  in  a  cubic  patch 
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Figure  6:  Avoiding  extra  sheets  in  a  cubic  patch 

5.  More  dramatically,  there  may  be  even  handles  etc.  on  a  cubic  patch.  See  Figure  3. 

Notice  that  we  listed  singular  points  together  with  self  intersection  because  for  implicit  patches, 
singular  points  appear  where  self  intersections  occur. 

We  use  tetrahedra  [X1X2WX4]  in  Figure  2  as  an  example  to  explain  how  the  situations  listed 
above  can  be  avoided  by  controlling  the  free  control  points  we  have  chosen.  I:.  this  example,  the 
free  control  points  are  the  function  value  h.i(x4)  and  the  gradient  'Vh^(x4).  The  same  argument 
with  minor  modifications  applies  to  the  cubic  polynomials  defined  for  other  tetrahedra. 

Situation  one  can  be  avoided  by  properly  choosing  the  function  value  at  x^.  Consider  the  line 
segment  from  the  centroid  of  [X1X2W],  p,  to  x.|.  If  the  function  value  fi3(x4)  is  chosen  to  be  have 
a  sign  opposite  to  that  of  the  function  value  at  p,  then  there  must  be  a  point  on  the  line  segment 
[x<jp]  where  the  cubic  polynomial  is  zero.  In  other  words,  the  zero  contour  passes  through  the 
interior  of  the  tetrahedron  [x^wx-jj. 

Situations  two  through  five  can  be  avoided  by  enforcing  monotonicity  conditions  on  the  cubic 
polynomial  along  the  direction  from  w  to  x,t.  A  function  is  monotone  in  direction  a  if  the  directional 
derivative  along  a  is  positive.  Let  the  cubic  polynomial  in  [xiX2X,|w]  be 

M"3)  =  Y,  «**’(**)• 

|A|=3 

A  sufficient  condition  for  the  cubic  polynomial  to  be  monotone  along  the  direction  form  w  to 
x.)  within  the  tetrahedron  is  that 

«A_ei+c<  -  «a  >  0.  for  all  A  with  Aj  >  1.  (18) 
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•  Location  of  shape  parameter 
Figure  7:  The  shape  control  scheme 

When  Aj  >  1,  the  condition  (18)  can  be  enforced  by  the  function  value  and  gradient  at  X4.  As  for 
Ai  =  1,  the  control  points  involved  in  (18)  are  completely  determined  from  the  prescribed  normals 
in  the  input  data,  so  the  monotonicity  conditions  may  not  be  satisfied  for  certain  input  data  no 
matter  how  the  free  control  points  h3(x4)  and  V/i,j(x^)  are  chosen.  But  remember  the  prescribing 
normals  in  the  input  data  is  only  a  measure  to  control  the  behavior  of  each  cubic  patch.  If  we 
choose  these  normals  within  proper  ranges,  the  condition  (18)  can  be  enforced. 

In  practice,  the  free  control  points  are  computed  using  the  degree- elevation  property  of  Bernstein- 
Bezier  representation.  The  idea  is  to  extend  the  effects  of  prescribed  normals  to  the  free  control 
points.  A  quadric  polynomial  q  over  the  tetiahedron  [X1X2X3X4]  can  be  determined  from  the  fact 

q{xi)  =  0,  V<7(x.)  =  n,-,  i  =  1,2.3. 

and  the  value  q(x4)  which  is  referred  to  as  a  shape  parameter.  If  this  is  done  for  all  facets  of 
the  input  polyhedron  V ,  then  quadric  polynomials  over  tetrahedra  such  as  [x2X3x"x.j]  can  be 
determined  also.  These  quadric  polynomials  aie  then  degree  elevated  to  cubic  polynomials,  whose 
control  points  corresponding  to  the  free  contiol  points  are  given  to  the  free  control  points.  This 
method  of  choosing  free  control  points  works  very  well  in  practice.  From  our  experience,  the  ranges 
of  free  control  points  within  which  the  cubic  patches  behave  well  are  fairly  large.  As  long  as  the 
free  control  points  are  not  in  the  relative  small  "bad”  ranges,  the  cubic  patches  are  in  good  shape. 

Figure  5  and  Figure  6  are  two  examples  of  how  the  above  method  works  in  the  setting  of  Figure 
2.  In  Figure  6,  the  left  figure  has  an  extra  sheet  due  to  badly  chosen  free  control  points.  In  the 
right  figure,  the  badly  chosen  free  control  points  are  corrected  using  the  above  method.  Figure  5  is 
similar  except  the  problem  is  the  hole  in  the  left  figure. 

3.4  Features  of  the  algorithm 

The  above  free-form  algorithm  has  several  features.  From  the  description  of  the  algorithm,  it  is 
easy  to  see  that  the  free-form  construction  in  the  algorithm  is  local.  In  the  following,  we  discuss  the 
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•  Location  of  shape  parameter 


Figure  8:  Figure  7  with  a  shape  parameter  decreased 


quadratic  precision  of  the  algorithm,  and  how  to  control  the  shape  of  the  free  form  surface  without 
knowing  the  details  of  the  algorithm. 

Quadratic  precision  is  a  measure  of  accuracy  of  free  form  surface  algorithms  in  terms  of  how 
well  the  algorithms  can  reproduce  a  known  surface  if  the  input  data  is  taken  from  the  surface.  A 
question  that  users  often  ask  about  a  free  form  surface  algorithm  is  that  if  the  input  data  is  taken 
from  a  sphere,  can  the  algorithm  reproduces  the  sphere.  For  the  algorithm  we  derive,  the  answer 
is  yes.  In  fact,  the  algorithm  reproduces  all  quadrics. 

Notice  that  the  input  polyhedron  V,  prescribed  normals  at  the  vertices  of  V,  and  the  shape 
parameters  completely  determine  the  free  form  surface.  If  the  input  data  is  taken  from  a  quadric 
surface  and  the  shape  parameters  are  frortr  the  quadric  surface,  then  the  algorithm  will  produce 
the  same  quadric  surface.  To  ensure  the  shape  parameters  are  properly  chosen  so  that  all  quadric 
surfaces  can  be  reproduced,  we  must  give  certain  default  values  to  the  shape  parameters.  For 
example,  an  easy  way  to  do  so  is  as  follows.  Randomly  choose  enough  vertices  of  V  so  that  these 
vertices  determine  a  quadratic  polynomial  q  such  that  the  zero  contour  of  q  passes  though  the 
chosen  vertices,  then  compute  the  the  shape  parameters  by  evaluating  q. 

An  important  feature  of  the  algorithm  we  derive  is  that  it  allows  the  users  to  control  the  shapes 
of  the  free-form  surfaces  produced  by  the  algorithm  without  knowing  the  details  of  the  algorithm. 
This  feature  is  very  important  for  applications  like  CAD/CAM,  where  the  designers  manipulate 
the  shape  of  the  free-form  surfaces  to  achieve  functional  or  aesthetic  design  objectives. 

Recall  that  for  each  facet,  the  cubic  polynomials  over  the  double  tetrahedron  containing  the 
facet  is  not  completely  fixed.  A  double  tetrahedron  has  a  vertex  outside  V,  and  we  call  the  vertex 
the  apex  of  the  double  tetrahedron.  At  the  apex  of  the  double  tetrahedron  of  each  facet,  the  value 
of  the  cubic  polynomials  is  left  as  a  shape  parameter,  as  was  shown  in  3.3. 

If  we  think  of  the  algorithm  as  producing  the  global  function  /  over  the  constructed  neighbor¬ 
hood  S  of  V,  then  the  value  of  /  at  each  apex  is  a  shape  parameter.  Since  the  interior  of  the 
free-form  surface  is  exactly  the  region  where  /  <  0.  decreasing  a  shape  parameter  at  a  apex  pulls 
the  free-form  surface  towards  the  apex.  Moreover,  only  nearby  quadric  patches  are  affected  by 


Figure  9:  An  example  of  shape  control 

this  shape  parameter  because  the  free-form  surface  construction  in  the  algorithm  is  local.  So,  the 
apexes  form  a  net  which  controls  the  shape  of  the  free-form  surface  through  the  she.-  -  parameters 
at  the  apexes. 

Figure  7  and  Figure  8  illustrate  a  two  dimensional  analogy  of  this  shape  control  scheme.  The 
situation  in  the  three  dimension  is  the  same  but  harder  to  draw.  Figure  9  is  an  example  of  two 
free-form  surface  having  everything  identical  except  the  shape  parameters. 

4  Conclusions 

We  have  presented  an  algorithm  for  generating  fiec-form  surfaces  from  the  input  data  of  an  arbitrary 
polyhedron  with  a  normal  prescribed  at  each  veitex  of  the  polyhedron.  The  algorithm  constructs  a 
smooth  piecewise  cubic  surface  interpolating  the  vertices  of  the  input  polyhedron  and  the  prescribed 
normal  at  each  vertex.  The  free-form  surface  construction  is  local  and  quadratically  precise.  In 
addition,  the  free-form  surface  produced  can  be  manipulated  through  a  set  of  intuitive  shape 
parameters  without  knowing  the  details  of  the  algorithm. 
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Figure  10:  A  skewed  dodecahedron 


1C 


Figure  11:  A  tea  pot 
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Figure  10  and  Figure  11  illustrate  some  implemented  results.  Figure  10  is  a  skewed  dodecahe¬ 
dron  with  12  points,  20  facets,  and  80  patches;  Figure  11  is  a  tea  pot  with  45  points,  72  facets,  and 
266  patches.  These  two  pictures,  as  well  as  the  pictures  shown  earlier,  are  generated  by  polygonizing 
the  cubic  patches  and  rendering  the  resultant  polygon  using  Gouraud  shading. 

We  hope  to  incorporate  the  free-form  surface  algorithm  into  a  geometric  modeling  system  and 
to  experiment  designing,  manipulating,  and  reasoning  about  complex  smooth  objects. 
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Abstract 

The  recent  explosion  of  interest  in  physical  system  simulation  may 
soon  lead  to  realistic  animation  of  passive  objects,  such  as  sliding  blocks 
or  bouncing  balls.  However,  complex  active  objects  (like  human  figures 
and  insects)  need  a  control  mechanism  to  direct  their  movements.  We 
present  a  paradigm  that  combines  the  advantages  of  physical  simula¬ 
tion  and  algorithmic  specification  of  movement.  The  animator  writes 
an  algorithm  to  control  the  object  and  runs  this  algorithm  on  a  phys¬ 
ical  simulator  to  produce  the  animation.  Algorithms  can  be  reused  or 
combined  to  produce  complex  sequences  of  movements,  eliminating  the 
need  for  tedious  keyframing.  We  have  applied  this  paradigm  to  control 
a  walking  biped.  The  walking  algorithm  is  presented  along  with  the 
results  from  testing  with  the  Newton  simulation  system. 


1  Introduction 

This  paper  describes  a  new  paradigm  for  the  control  and  animation  of  complex 
active  objects  such  as  the  human  figure.  This  opproach  allows  the  animator 
to  control  an  object  through  an  algorithm  which  specifies  certain  “intuitive” 
variables  as  a  function  of  time  and  of  world  state.  In  the  case  of  human  figure 
walking,  the  animator  might  write  an  algorithm  which  controls  the  acceleration 
of  the  figure’s  center  of  mass  at  one  point  in  the  animation,  and  which  controls 
the  angle  of  the  knees  at  another  point.  The  algorithmic  approach  to  animation 
allows  this  to  be  done  with  ease,  as  demonstrated  by  the  walking  algorithm 
presented  in  Section  6. 
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Witkin  and  Kass  (WK88|  have  combined  physical  simulation  and  key¬ 
framing  to  produce  realistic  animation  of  their  jumping  Luxo  lamp.  With 
their  approach  the  animator  uses  spacetime  constraints  to  specify  several  key 
points  for  selected  variables.  These  variables  may  be  positions,  velocities, 
forces  and  so  on.  Combining  spacetime  constraint  equations  with  the  La- 
grangian  equations  of  motion  and  discretizing  over  time  yields  a  system  of 
equations  that  are  solved  to  produce  the  motion.  Since  the  system  is  gen¬ 
erally  underconstrained  (having  multiple  solutions)  a  solution  can  be  chosen 
to  minimize  the  power,  fuel  comsumption  and  so  on. 

Our  algorithmic  approach  is  similar  in  that  the  animator  can  control 
accelerations  and  forces,  but  differs  in  that  the  constraints  can  be  added 
or  removed  “on  the  fly”  as  the  algorithm  sees  changes  in  the  world  state 
which  might  not  be  predictable.  In  the  case  of  human  figure  walking  the 
algorithm  might,  as  the  foot  touches  the  ground,  remove  a  foot  positioning 
constraint  and  add  a  leg  stiffening  constraint.  The  exact  point  of  contact 
is  not  predictable  in  advance.  Additionally,  the  algorithmic  approach  frees 
the  animator  from  considering  the  dynamics  of  impact  and  other  changes  in 
kinematic  relationships,  which  are  handled  automatically  by  the  simulation 
component  of  our  system.  Incorporating  impact  into  the  work  of  Witkin 
and  Kass  would  require  either  guessing  the  impact  points  beforehand  or 
incorporating  a  “force  field”  approach  as  described  in  Section  2. 

Other  work  on  combining  control  and  simulation  has  been  done  by  Barzel 
and  Barr  (BB88).  Their  method  of  dynamic  constraints  adds  fictitious  forces 
which  pull  the  simulated  objects  into  specified  positions.  By  doing  this  in  the 
framework  of  a  simulation  system,  the  movement  of  complex  physical  objects 
can  be  simulated  with  little  work  on  the  part  of  the  animator.  A  limited  form 
of  control  is  achieved  by  attaching  forces  to  points  on  the  object  and  dragging 
these  points. 

Various  other  approaches  to  combine  control  and  physical  simulation  have 
been  explored.  Wilhelms  (Wil87j  blends  kinematic  and  dynamic  formula¬ 
tions,  Isaacs  and  Cohen  [IC87]  incorporate  inverse  dynamics  in  their  simula¬ 
tion  system,  and  Brotman  and  N' travali  [BN881  use  dynamics  and  optimal 
control  to  interpolate  between  key  frames. 

Some  further  ins.ghts  on  control  can  be  gained  from  examining  the  current 
literature  in  the  field  of  robotics.  While  this  field  deals  with  controlling  real, 
physical  objects,  some  of  the  techniques  can  be  applied  to  produce  simpler 
animation. 

Researchers  in  robotics  have  taken  various  approaches  to  reduce  the  com¬ 
plexity  of  control  programs  for  physical  objects.  The  computed  torque 
method  (see  (Cra86|)  for  robot  arms  can  be  viewed  as  simplifying  control 
by  reducing  the  gripper  to  a  unit  mass.  The  control  program  can  ignore  the 
dynamics  of  the  robot  arm,  only  concerning  itself  with  the  position  of  the 
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end  effector  as  a  function  of  time. 

In  building  his  one-legged  hopping  machine,  Raibert  [Rai86j  partitioned 
control  along  three  intuitive  degrees  of  freedom:  hopping,  forward  speed  and 
body  posture.  This  resulted  in  surprisingly  simple  control  programs  for  the 
hopping  robot.  For  multi-legged  machines,  Raibert  introduced  the  idea  of  a 
“virtual  leg”  which  was  defined  in  terms  of  the  robot’s  physical  legs.  This 
again  led  to  simplified  control  programs. 

Both  the  computed  torque  method  and  Raibert’s  virtual  leg  demonstrate 
that  a  proper  choice  of  control  variables  can  lead  to  simplified  control  pro¬ 
grams.  The  problem  with  this  approach  is  that  there  is  often  no  simple 
closed-form  mapping  of  these  control  variables  omo  the  forces  and  torques 
needed  to  control  the  object.  In  some  cases  a  complete  system  of  equations 
must  be  numerically  solved  to  make  this  mapping.  This  is  called  “inverse 
dynamics”  and  is  typically  rejected  by  robotics  researchers  as  being  too  ex¬ 
pensive  to  use  in  real-time  control.  For  the  purposes  of  animation,  however, 
it  is  ideal. 

This  is  the  basis  of  our  algorithmic  approach  to  control.  This  approach 
advocates  the  selection  of  a  small  set  of  intuitive  variables  which  are  used 
by  the  algorithm  in  controlling  the  object.  The  algorithm  constrains  these 
variable  with  constraint  equations ,  which,  when  combined  with  the  standard 
Newton-Euler  equations  of  motion,  produce  a  system  of  equations  describing 
the  motion  of  the  simulated  object.  The  system  of  equations  is  maintained 
by  our  general  purpose  physical  simulator,  called  Newton.  The  Newton  sim¬ 
ulator  is  responsible  for  integrating  the  motion  of  the  simulated  objects  over 
time  to  produce  the  animation.  As  described  in  the  next  section,  Newton 
also  automatically  updates  the  system  of  equations  as  kinematic  relation¬ 
ships  in  the  simulation  change  (one  such  change  would  occur  as  the  biped’s 
foot  touches  the  ground).  Finally,  Newton  provides  an  interface  to  allow  the 
algorithm  to  add  and  remove  constraint  equations  to  and  from  the  system 
of  motion  equations. 

In  the  event  that  the  control  algorithm  underconstrains  the  motion  of 
the  object,  constrained  optimization  techniques  are  used  to  choose  a  motion 
that  optimizes  some  criierion  while  satisfying  the  constraints  imposed  by  the 
algorithm.  Our  decision  to  allow  control  programs  to  underconstrain  the  con¬ 
trolled  object  -  necessitating  the  use  of  constrained  optimization  techniques 
-  is  based  on  the  realization  that  control  algorithms  often  require  many  fewer 
control  variables  than  there  are  degrees  of  freedom  in  the  controlled  object. 
A  robot  modeled  after  the  human  figure  may  have  as  many  as  two  hundred 
degrees  of  freedom  [Zel82j,  while  the  control  program  for  such  a  robot  would 
only  require  twenty  or  thirty  degrees  of  freedom  to  accomplish  its  task.  In 
programming  our  walking  biped  we  used  at  most  eleven  of  its  sixteen  degrees 
of  freedom  at  any  given  instant. 
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In  summary,  the  algorithmic  approach  presented  in  this  paper  allows  the 
algorithm  to  constrain  a  small  set  of  intuitive  variables.  The  algorithm  is 
allowed  to  underconstrain  the  motion  of  the  object,  in  which  case  a  motion 
is  chosen  which  optimizes  some  criterion  while  obeying  the  constraints.  The 
Newton  simulator  incorporates  the  constraint  equations  into  its  automat¬ 
ically  maintained  system  of  motion  equations  and  integrates  over  time  to 
produce  realistic  animation. 

Section  2  outlines  the  relevant  background  of  the  Newton  simulation  sys¬ 
tem.  Section  3  describes  in  detail  the  algorithmic  approach,  while  Section  4 
looks  at  some  low-level  controllers  used  by  the  walking  algorithm.  Following 
this,  Sections  5  and  6  outline  the  biped  model  and  the  walking  algorithm, 
and  present  results  from  testing  the  algorithm. 

2  Overview  of  Newton 

The  walking  algorithm  described  in  this  paper  has  been  designed  and  tested 
using  the  Newton  simulation  system,  part  of  a  large  research  effort  in  mod¬ 
eling  and  simulation  at  Cornell  University.  The  development  of  Newton  was 
inspired  by  the  need  for  more  general-purpose,  flexible  simulation  systems. 

Extensive  mechanical  engineering  research  has  led  to  many  developments 
in  physical  system  simulation.  The  ADAMS  (Cha.85)  and  DADS  (HL87| 
systems  are  examples  of  large  state-of-the-art  systems  from  the  mechani¬ 
cal  engineering  domain.  In  many  ways  such  systems  are  very  sophisticated: 
efficient  formulations  of  mechanism  dynamics  are  supported,  fancy  numer¬ 
ical  techniques  for  solving  equation  systems  are  used,  object  flexibility  and 
elasticity  are  often  handled,  and  so  on.  Recent  work  by  graphics  and  ani¬ 
mation  researchers  [BB88,IC87,MW88,Hah88]  in  what  is  termed  physically- 
based  modeling  has  generally  been  less  sophisticated  but  has  placed  greater 
emphasis  on  animation  of  interesting  high-degree-of-freedom  mechanisms. 

A  number  of  things  are  still  lacking  in  all  of  these  systems.  Typically  they 
have  almost  ignored  geometric  considerations  and  represented  objects  simply 
as  point  masses  with  associated  inertias  and  coordinate  systems.  Geometric 
modeling  techniques  have  matured  enough  to  allow  object  representations 
used  by  dynamic  simulations  to  include  a  complete  geometric  description 
usable  by  a  geometry  processing  module.  Furthermore,  impact,  contact,  and 
friction  are  typically  handled  by  current  systems  in  an  ad  hoc  or  rudimentary 
manner,  if  at  all.  In  some  cases,  for  instance,  any  possible  impacts  must  be 
specified  in  advance:  in  others,  a  kind  of  “force  field”  technique  is  used,  in 
which  between  every  pair  of  objects  there  is  a  repelling  force  that  is  negligible 
except  when  objects  are  very  close  together.  In  addition,  the  desire  to  manip¬ 
ulate  high-degree-of-freedom  objects  suggests  that  a  module  for  specification 
of  control  algorithms  should  be  a  significant  part  of  a  dynamics  system. 


4 


2.1  Newton  Architecture 


Using  Newton,  a  designer  can  define  complex  three-dimensional  physical  ob¬ 
jects  and  mechanisms  and  can  represent  object  characteristics  from  a  wide 
range  of  domains.  An  object  is  made  up  of  a  number  of  ‘‘models,’’  each 
responsible  for  organization  of  object  characteristics  from  a  particular  do¬ 
main.  In  most  simulations  the  basic  domains  of  geometry,  dynamics,  and 
controlled  behavior  are  modeled.  A  dynamic  modeling  system,  for  example, 
is  responsible  for  maintaining  an  object’s  position,  velocity,  and  accelera¬ 
tion,  and  for  automatically  formulating  the  object’s  dynamics  equations  of 
motion.  A  geometric  modeling  system  is  responsible  for  information  about 
an  object’s  shape,  distinguished  features  on  the  object,  and  computation 
of  geometric  integral  properties  such  as  volume  and  moments  of  inertia.  It 
also  detects  and  analyzes  obje-t  interpenetrations  so  that  an  interference 
modeling  system  can  deal  with  collisions  between  objects. 

Newton  is  composed  of  three  main  components:  the  definition  and  repre¬ 
sentation  module,  the  analysis  module  and  the  report  system.  The  definition 
module  analyzes  high  level  language  descriptions  of  Newton  entities  and  orga¬ 
nizes  the  corresponding  data  structures.  The  analysis  component  implements 
the  top-level  control  loop  of  simulations  and  coordinates  the  working  of  vari¬ 
ous  analysis  subsystems.  The  report  system  handles  generation  of  graphical 
feedback  to  users  during  simulations  as  well  a s  recording  of  relevant  infor¬ 
mation  for  later  regeneration  of  animations. 

2.2  Dynamic  Analysis  in  Newton 

A  comp'  -x  physical  object  is  modeled  as  a  collection  of  rigid  bodies  related 
by  constraints.  Newton- Euler  equations  of  motion  are  associated  with  each 
individual  rigid  body.1  At  the  time  an  object  is  created  the  equations  are  of 
the  form 

mf  =  0 

Jw  -r  u;  x  Jw  =  0. 

where  m  is  the  mass,  r  is  the  second  time  derivative  of  the  position  (ie.  the 
acceleration),  J  is  the  3x3  inertia  matrix,  and  w  and  u;  are  the  rotational 
velocity  and  acceleration,  respectively. 

A  specification  that  two  objects  are  to  be  connected  with  a  spherical 
hinge  is  met  by  the  addition  of  one  vectorial  constraint  equation  and  the 
addition  of  some  terms  to  the  motion  equations  of  the  constrained  objects. 
For  a  holonomic  constraint  such  as  this  one,  the  second  derivative  of  the 
constraint  equation  can  be  used  along  with  the  modified  motion  equations 

1  Newton  is  capable  of  using  dynamics  formulations  other  than  the  one  outlined  here.  We 
are  also  working  on  incorporating  non-rigid  bodies  into  the  system. 


to  solve  for  object  accelerations  and  reaction  forces.  Thus,  the  equations 
above  become 

—  P hinge 

J  1^1  rta)|  A  J[u>i  =  C[  X  Fh,ngc 

TTI2T2  —  F kt  ngc 

2  T  u72  X  J2 u72  =  C2  X  -Fhmgc 

7*1  4"  Uq  X  C(  4-  U>i  X  (uq  X  C[  )  =  7*2  4”  u72  X  C2  4"  uJ2  X  (u72  X  C2), 

where  c,  is  the  vector  from  object  i’s  center  of  mass  to  the  location  of  the 
hinge  and  Fh,ngc  is  the  constraint  force  that  keeps  the  objects  together.  Note 
that  the  last  equation  above  is  the  second  time  derivative  of  the  holonotnic 
constraint  equation  rq  4-  ct  =  r2  -r  c2  for  spherical  joints.  Other  kinds  of 
hinges  commonly  used  in  Newton  include  revolute  or  pin  joints,  prismatic 
joints,  springs  and  dampers,  and  rolling  contacts. 

If  gravity  is  present  during  the  simulation  the  system  will  automatically 
add  gravitational  force  terms  to  the  objects’  translational  motion  equations. 
The  system  keeps  track  of  the  constraints  responsible  for  the  various  terms 
in  the  motion  equations.  Thus,  constraints,  and  their  corresponding  motion 
equation  terms,  can  be  removed  at  any  time  without  necessitating  complete 
rederivation  of  the  system  of  motion  equations. 

Using  this  method  of  dynamics  formulation,  closed-loop  kinematic  chains 
are  handled  as  simply  as  open  chains.  Though  the  formulation  does  lead  to 
a  large  set  of  equations,  the  matrices  are  very  sparse  and  often  symmetric. 
Thus,  acceptable  efficiency  is  achieved  by  the  use  of  sparse  matrix  solution 
techniques. 


2.3  Event  handling,  impact  and  contact 

Newton,  unlike  many  other  simulation  systems  (though  see  [Fea85;),  can 
automatically  and  incrementally  reformulate  the  motion  equations  as  excep¬ 
tional  events  occur  during  simulations.  One  kind  of  exceptional  event  is  a 
change  in  kinematic  relationship  between  objects.  Figure  1  shows  a  block 
that  was  initially  sliding  along  a  table  top.  After  some  time  the  edge  of 
the  table  is  reached  and  the  contact  relationship  changes  from  a  plane-plane 
contact  to  a  plane-edge  contact.  Still  later  the  contact  is  broken  altogether. 
These  changing  contact  relationships  are  automatically  detected  by  Newton. 
The  system  of  motion  equations  and  the  related  constraint  equations  are 
automatically  maintained  by  Newton  to  reflect  these  changing  relationships. 

During  the  course  of  a  simulation,  a  variety  of  events  can  occur  that 
require  special  processing.  Newtons  event  handler  is  primarily  responsible 
for  detection  and  resolution  of  impacts,  for  analysis  of  continuous  contacts 
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Figure  1:  Changing  Kinematic  Relationships 

between  objects  and  corresponding  maintenance  of  temporary  hinges ,  special 
kinds  of  hinges  that  model  one  sided  constraints  between  objects  in  contact, 
and  for  handling  of  events  specified  by  control  programs  that  necessitate 
changes  in  the  constraint  set.  For  example,  the  walking  algorithm  might  tell 
the  event  handler  to  notify  it  when  the  biped’s  foot  touches  the  ground  so 
that  it  can  change  the  constraint  equations. 

The  geometric  modeling  subsystem  is  responsible  for  detecting  and  an¬ 
alyzing  impacts  and  interpenetrations.  In  the  usual  method  of  handling 
impacts,  the  dynamic  analysis  module  formulates  impulse-momentum  equa¬ 
tions  in  a  manner  completely  analagous  to  the  formulation  of  the  basic 
dynamics  equations,  and  solves  these  equations  to  produce  the  instanta¬ 
neous  velocity  changes  caused  by  the  impact.  The  details  of  Newton's  meth¬ 
ods  for  handling  impact,  contact  and  other  exceptional  events  are  given  in 
[HH87,HH88,CS88,Cre89]. 

3  The  Algorithmic  Approach 

In  Newton’s  automatically-generated  equations  of  motion  certain  quantities 
are  considered  to  be  unknovms.  A  system  of  simultaneous  linear  equations  is 
solved  at  each  time  step  to  produce  values  for  the  unknowns.  These  values 
are  integrated  over  time  to  produce  the  simulated  motion.  Typically,  the 
unknowns  consist  of  accelerations  and  joint  constraint  forces,  while  positions, 
velocities  and  joint  control  torques  are  hnoums. 

In  the  algorithmic  approach,  the  programmer  controls  “intuitive”  quanti¬ 
ties  defined  as  linear  combinations  of  the  unknowns.  The  programmer  might, 
for  example,  want  to  control  the  acceleration  of  the  center  of  mass  of  a  biped 
without  explicitly  controlling  each  component  of  the  biped.  To  do  this,  the 
algorithm  must  define  the  acceleration  of  the  center  of  mass  in  terms  of  the 
accelerations  of  the  centers  of  mass  of  the  primitive  components  of  the  ob- 


procedure  initialize 
begin 

add-equation  "  fcm  =  ^  ^  m,  f,  " 
end 

procedure  controller(  time  ) 
begin 

=  /{  time  ) 
end 


Figure  2:  The  Format  of  an  Algorithm 

ject.  Over  the  course  of  execution,  the  algorithm  must  supply  the  desired 
acceleration  of  the  center  of  mass  at  each  point  in  time. 

Figure  2  shows  the  format  of  a  control  algorithm.  For  the  sake  of  clarity 
the  algorithms  will  be  described  in  a  Pascal-like  notation2.  Two  procedures 
are  always  present:  one  to  initialize  the  algorithm  (called  initialize)  and 
one  to  be  executed  repeatedly  over  the  course  of  the  task  (called  controller). 
The  controller  procedure  has  access  to  the  complete  state  of  the  system. 
The  algorithm  of  Figure  2  trivially  defines  and  controls  the  acceleration  of 
the  center  of  mass  of  an  object  (the  function  /  must  be  defined  elsewhere). 

Defining  and  controlling  a  three-dimensional  vectorial  quantity  like  the 
acceleration  of  the  center  of  mass  has  the  effect  of  adding  three  constraint 
equations  to  the  system  of  simultaneous  Unear  equations  that  describe  the  in¬ 
stantaneous  motion  of  the  object.  By  considering  joint  torques  as  unknowns 
in  this  augmented  system  of  equations,  the  system  can  be  solved  to  produce 
motion  that  satisfies  the  additional  constraint  equations.  This  is  a  simple 
apphcation  of  inverse  dynamics. 

For  an  object  with  n  degrees  of  freedom  the  control  algorithm  can  define 
and  control  up  to  n  independent  scalar  quantities3.  If  fewer  than  n  equations 
are  added  the  system  of  motion  equations  is  underdetermined,  and  many  dif¬ 
ferent  solutions  could  satisfy  the  constraints  of  the  control  algorithm.  In  this 
case  the  algorithm  must  guide  the  selection  of  a  solution  by  providing  a 
cost  function  which  is  quadratic  in  the  unknowns.  A  standard  numerical 
optimization  technique  is  used  to  compute  a  solution  that  instantaneously 
(for  each  point  in  time)  minimizes  the  cost  function  while  obeying  the  algo¬ 
rithm’s  constraints.  This  is  different  from  the  approach  of  Witkin  and  Kass 

2The  algorithms  are,  for  now,  written  in  Lisp. 

3The  additional  definitional  equations  could  make  the  system  of  motion  equations  incon¬ 
sistent.  This  would  be  an  error  on  the  part  of  the  control  algorithm. 
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[WK88j,  who  optimize  over  the  whole  animation.  This  reflects  the  different 
philosophies  of  the  two  systems:  Witkin  and  Kass  specify  all  of  the  infor¬ 
mation  beforehand,  while  we  let  the  control  algorithm  make  decisions  during 
the  animation.  Such  “on  the  fly”  decisions  make  it  impossible  to  do  global 
optimization,  but  allow  much  more  versatility  in  the  control  algorithm  by 
not  requiring  a  priori  knowledge  of  impacts  and  other  exceptional  events. 

In  summary,  the  programmer  designs  an  algorithm  in  a  high-level  com¬ 
puter  language  to  control  intuitive  degrees  of  freedom  of  the  object.  These 
degrees  of  freedom  are  defined  as  linear  combinations  of  the  unknowns  in 
the  object’s  equations  of  motion.  An  augmented  linear  system  of  equations 
describes  the  instantaneous  behavior  of  the  object;  this  system  can  be  solved 
to  produce  the  object’s  configuration  at  each  point  in  time.  If  the  system 
is  underdetermined,  the  algorithm  can  provide  a  cost  function  to  guide  the 
choice  of  a  solution. 

In  the  remaining  sections  we  describe  the  application  of  this  approach  to 
the  design  of  a  simple  walking  algorithm. 


4  Low-level  Controllers 

In  designing  algorithms  with  Newton  we  found  ourselves  frequently  using  PD 
controllers'*  and  curve-fitting  controllers  to  control  the  “trajectory”  of  many 
of  the  defined  quantities.  In  controlling  the  biped,  for  example,  a  quintic 
interpolation  was  used  to  plot  the  trajectory  of  the  heel,  and  a  PD  controller 
was  used  to  orient  the  foot  before  it  struck  the  ground.  A  small  library  of 
these  controllers  is  used  in  the  biped  algorithm,  and  will  be  described  here. 

PD  controllers  are  used  in  the  biped  algorithm  to  control  orientation, 
position  and  joint  angle.  Each  controller  adds  an  equation  to  the  system 
of  motion  equations  which  defines  the  second  derivative  of  the  quantity  in 
terms  of  the  first  derivative  and  the  quantity  itself.  The  procedure  in  Fig¬ 
ure  3  produces  accelerations  to  move  an  object  to  within  1%  of  a  position 
x-desired  within  a  given  time  delta-time.  The  quantities  x,  v  and  a  are 
data  structures  representing  state  variables  of  the  controlled  object.  These 
data  structures  are  used  by  the  add-named- equation  function  to  create  the 
appropriate  equation. 

Execution  of  the  procedure  in  Figure  3  causes  a  named  equation  to  be 

*A  PD  controller  (Proportional,  Derivative),  also  known  as  a  “spring  and  damper”  con¬ 
troller,  relates  the  second  derivative  of  a  variable  linearly  to  the  error  in  the  variable’s  first 
derivative  and  to  the  error  in  the  variable  itself.  The  equation  is  £+  -  Zd«ir<d)  =  0 

for  some  appropriate  r.  PD  controllers  are  used  extensively  in  robotics  to  move  robot  joints 
into  specified  positions  by  calculating  the  joint  acceleration  as  a  function  of  the  position  and 
velocity  errors.  A  good  explanation  can  by  found  in  [Cra86|.  Barzel  and  Barr  [BB88]  use  a 
form  of  PD  controller  to  achieve  their  dynamic  constraints. 
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procedure  position-«ith-PD(  constraint-name ,  object, 

x-desired,  delta-time  ) 


var  z,  v,  a:  quantity 
t:  real 

begin 

z  =  get-position-quantity(  object  ) 
v  =  get-velocity-quantity (  object  ) 
a  =  get-acceleration-quantity(  object  ) 

t  =  -  delta-time  /  log(  .01  ) 

add-naaed-equation(  constraint-name , 

"  a  r  ;  u-r  ;r( z  -  x-desired)  =  0  ”  ) 

end 


Figure  3:  PD  Controller  Used  in  Positioning 

added  to  the  system  of  motion  equations.  This  equation  will  continue  to 
affect  the  motion  of  the  object  until  it  is  explicitly  removed  by  the  control 
algorithm. 

A  complete  list  of  controllers  available  to  the  biped  walking  algorithm 
is  shown  in  Figure  7  at  the  end  of  the  paper.  Those  with  quintic  in  their 
name  do  quintic  interpolation  to  achieve  the  desired  position  and  velocity  in 
the  desired  time.  Quintic  interpolation  was  chosen  over  cubic  interpolation 
to  eliminate  “jerk”  (discontinuous  acceleration)  from  the  beginning  and  end 
of  the  trajectory. 


5  The  Biped  Model 

The  simulated  biped  is  composed  of  a  torso,  two  legs  with  knee  joints  and  two 
feet  with  toe  joints.  This  model,  was  adapted  from  a  description  in  (McM84] 
and  is  shewn  in  Figure  4.  The  hips  and  ankles  are  three  degree  of  freedom 
spherical  joints,  while  the  knees  and  toes  are  one  degree  of  freedom  revolute 
joints,  making  a  total  of  sixteen  degrees  of  freedom.  The  biped  is  about  six 
feet  tall  with  moments  approximating  those  of  a  human  being. 

We  hope  to  improve  this  model  by  incorporating  joint  limits  and  elas¬ 
tic  tendons.  McMahon  suggests  that,  during  walking,  energy  is  stored  in 
stretched  tendons  and  is  released  when  the  stretched  leg  swings  forward 
(McM84j.  This  idea  might  be  used  to  simplify  the  walking  algorithm  de¬ 
scribed  in  the  next  section. 

Newton's  impact  handling  capabilities  have  not  yet  been  extended  to 
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accurately  model  the  impact  of  the  feet  upon  the  ground.  Instead,  impact  is 
simulated  by  adding  an  external  force  and  torque  to  the  feet  that  holds  them 
level  with  the  ground  until  they  are  released  with  an  explicit  command  from 
the  control  algorithm.  This  is  as  though  the  biped  was  walking  with  magnetic 
shoes  on  a  steel  plate.  Very  shortly  we  expect  to  adapt  the  algorithm  to 
incorporate  realistic  impact. 

6  The  Walking  Algorithm 

An  abbreviated  version  of  the  walking  algorithm  is  shown  in  Figures  8  and 
9,  which  can  be  found  at  the  end  of  this  paper.  The  algorithm  cycles 
through  a  set  of  six  states:  swing  the  right  leg,  land  the  right  foot,  lift 
the  left  foot,  swing  the  left  leg,  land  the  left  foot,  lift  the  right  foot  and 
then  repeat  the  cycle.  In  the  swing  phase,  a  quintic  trajectory  is  plot¬ 
ted  for  the  swing  foot  with  aove-heel-to-target,  while  the  stance  leg  is 
stiffened  with  set-angle-sith-PD  and  the  foot  is  oriented  for  landing  with 
orient-sith-PD  (shown  under  START  in  Figure  9).  In  the  landing  phase, 
the  leading  leg  is  stiffened  as  the  foot  nears  the  ground.  Following  this,  the 
takeoff  phase  flexes  the  trailing  leg,  causing  the  trailing  foot  to  lift  from  the 
ground.  Once  the  trailing  toe  is  bent  to  10°  the  flexing  constraint  is  removed 
and  the  swing  phase  begins  for  the  trailing  leg. 

The  largest  number  of  constraints  are  applied  during  the  swing  phase,  as 
shown  in  Table  1.  Since  the  biped  has  sixteen  degrees  of  freedom  (DOF)  it 
remains  underconstrained  at  all  times.  A  quadratic  cost  function  is  therefore 
defined  (in  initialize  of  Figure  9)  in  order  to  fully  determines  the  motion 
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Constraint  Name  DOF  Constrained  Item 


TORSO-CONSTRAINT  3 

L-KNEE-ANGLE  l 

R-HEEL-TRAJ  3 

R- FOOT-ORIENTATION  3 

R-TOE- ANGLE  t 


torso  orientation  in  3  dim 
angle  of  revolute  knee  joint 
heel  acceleration  in  3  dim 
foot  orientation  in  3  dim 
angle  of  revolute  toe  joint 


Table  1:  Swing  Phase  Constraints 


of  the  biped.  The  cost  function  is  a  weighted  sum  of  the  translational  and 
angular  accelerations,  and  of  the  difference  between  the  torso  translational 
acceleration  and  some  acceleration  defined  by  a  function  F  which  tries  to 
keep  the  torso  mid-way  between  the  two  feet. 

We  found  that  a  cost  function  which  minimizes  instantaneous  transla¬ 
tional  and  rotational  acceleration  usually  produces  smooth  motion.  In  the 
case  of  the  simulated  biped,  the  cost  function  causes  the  constrained  heel 
acceleration  to  be  achieved  by  a  linear  combination  of  small  accelerations  of 
many  components  of  the  body,  rather  than  a  few  large  accelerations  of  those 
components  which  are  near  the  heel.  We  have  observed  that  the  combina¬ 
tion  of  many  small  accelerations  yields  more  stable  motion  than  large,  local 
accelerations. 

The  walking  algorithm  was  tested  with  the  Newton  simulation  system. 
Figure  5  shows  ten  frames  in  which  the  biped  completes  a  full  cycle  of  the 
six  phases  described  above.  The  full  simulation  consisted  of  twenty  seconds 
of  straight-line  walking  on  a  fiat  surface  and  generated  the  statistics  shown 
in  Figure  6.  The  version  of  the  algorithm  that  produced  these  statistics  had 
the  biped  increase  speed  at  4.0  seconds,  as  can  be  seen  on  the  graphs. 
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7  Summary 

We  have  presented  an  algorithmic  approach  to  control.  This  approach  allows 
the  animator  to  choose  intuitive  degrees  of  freedom  by  which  to  control  an 
object.  The  control  algorithm  adds  and  removes  constraint  equations  “on 
the  fly”  as  the  world  state  changes;  a  priori  knowledge  of  the  exact  mo¬ 
ment  of  each  state  change  is  not  required.  With  the  algorithmic  approach, 
all  consideration  of  dynamics  and  impact  is  left  to  the  Newton  simulation 
system.  The  burden  on  the  animator  is  further  reduced  by  allowing  underde¬ 
termined  specification  of  motion  through  the  use  of  constrained  optimization 
techniques. 

We  have  presented  an  algorithm  to  control  a  simulated  biped,  along  with 
results  from  its  execution  on  the  Newton  simulation  system.  The  algorithm 
has  the  advantage  of  being  intuitive,  simple  to  program,  and  reusable. 

Unlike  keyframing,  the  algorithmic  approach  does  not  require  the  anima¬ 
tor  to  repeat  the  work  of  creating  new  key  frames  for  every  walking  sequence. 
Unlike  keyframing,  the  algorithmic  approach  allows  various  algorithms  to  be 
combined  to  produce  long  animated  sequences.  We  believe  that  in  the  future, 
animating  complex  physical  objects  will  require  a  structured,  algorithmic  ap¬ 
proach  similar  to  that  presented  in  this  paper. 

8  Future  Work 

We  will  incorporate  elastic  tendons  and  joint  friction  into  the  Newton  simu¬ 
lation  system  and  modify  the  walking  algorithm  accordingly.  From  there  we 
hope  to  develop  a  suite  of  algorithms  to  allow  a  biped  to  walk,  turn,  climb 
stairs,  manipulate  objects,  and  so  on.  In  keeping  with  the  structured  ap¬ 
proach  presented  in  this  paper  we  will  attempt  to  combine  these  algorithms 
to  have  the  biped  perform  complicated  tasks.  In  carrying  an  object  up  a 
flight  of  stairs  the  high-level  algorithm  would  combine  subroutines  to  pick 
up  the  object,  walk  to  the  stairs,  climb  the  stairs  and  deposit  the  object. 
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position-aith-PD(  constraint -nana,  objaet,  £d,  A;  5 

position-point  -»ith-PD(  constraint -nana ,  objaet,  poiuc-on-  objaet,  xd,  At  ) 
oriant-aith-PD(  constraint -nan*,  objaet,  4^,  At  ) 
sat-angla-tith-PD(  constraint-nan* ,  joint,  $d,  At  ) 

position-»ith-qumtic(  eonstraint-naaa ,  objaet,  xd,  vd,  At  ) 

position-pomt-Tith-qainvie(  eonstraint-au.,  objaet,  pomt-on-objact ,  xd,  vd,  At  ) 
oriant-gith-quintic(  constraint -nana ,  objact ,  <j>d,  4>  d  ,  At  ) 
sat-angla-with-qumtie( constraint-nama ,  joint,  Sd,  9d,  At  ) 


Figure  7:  Low-level  Controllers 


const 


tuaa-in-air 

strida 

diraction 

msida-stap-fraetion 

baal-Y-strika-spaad 

haal-I-strika-spaad 

foot-strilca-oriantation 

torso-onantation 


=  O.S  s 
=  O.S  ■ 

=(100) 

=  20  x 
=  -O.OS  »/s 
=  0.02  m/t 
-  10°  about  (001) 

*  -10”  about  (0  0  1) 


phasa:  itart  r-saing  r-land  1-lift  1-taiaof f  1-swing  X-land  r-lift  r-takaoff  ) 


procadnra  noTa-haal~to-targat(  constraint-nana ,  foot,  othar-foot,  hip,  othar-hip  ) 

»ar  targat-x,  targat-T,  hip-to-hip:  ractor 
bagin 

hip-to-hip  ■  gat-position(  TOESQ,  hip  )  -  gat-potition(  TOISO,  othar-hip  ) 

targat-x  =■  gat-position(  othar-foot,  BEJEL  )  ♦  strida  x  diraction 
♦  lasida-stap-fraction  x  hip-to-hip 

targat-T  =  haal-T-strike-  spaad  x  (0  t  0)  ♦  boa) -I-ctxlfcy-spcod  x  diraction 

posit  ion-point -vith-quint  ic  (  constraint -naao,  foot,  HTnt. ,  targat-x,  targat-T,  tma-in-a 
and 


Figure  8:  Definitions  for  the  Walking  Algorithm 


I 


procedure  initialize 

^  “  ^pl  rf-/aol)  "  P|a#i</I  ^  %-/ool)  "  ^(apjo) 

begin  _ 

quadratic-cost  -  p-  ^  r2  -  20  (ftvr<o  -  F)J 

phase  =  START 

end 

procedure  controlled  cut*  ) 
b«J« 

ease  phase  of 

START: 

phase  =  H-SWING 

orient-eith-PD(  TORSO-C-'f STRAIIT ,  TORSO,  torso-orientation,  2.0  s  ) 
aoee -heel- to- target!  g-BEEL-TRAJ .  R-BEEL,  L-BEEL,  R-HIP,  L-BIP  ) 
set-angle-sith-PD<  L-IIEE-AIGLE,  L-iJEE,  175°,  0.1  s  ) 

onent-»ith-PD(  R-FOOT-ORIEITATIOI ,  R-FGCT,  foot-strike-oneatation,  tiae-m-air  ) 
set-angle-«ith-PD(  R-TOE-ilGLE,  R-TOE-JOIIT,  0°,  tiae-m-air  ) 

R-SWING: 

if  distnnce-to-target(  R-FQOT  )  <  0.01  a  then 

phas.  =  R-LANDING 
renora-constraint (  R-BEEL-TRAJ  ) 

set-nngle-»ith-PD(  R-IIES-AIGLE,  R-XIEE,  t*S®,  O.OS  a  ) 

R-LANDING: 

if  heel-hes-touehed(  R-FOOT  )  than 
phas*  =  1,-TAKEOFF 

renoae-constraints!  R-FOOT-ORIEITATIOI.  R-TOE-AIGLE,  L-IIEE-AIGLE  ) 
set-angle-«ith-PD(  L-IIEE-AIGLE,  L-IIEE,  160°,  0.1  s  ) 

L-TAKEOFF: 

if  joint-angle!  L-TOE-JQIIT  )  >  10®  than 

phase  =  1,-SWING 
renora-constraint !  L-IIEE-AIGLE  ) 

aore-heel-to-target!  L-BEEL-TRAJ ,  L-IEXL,  I-IEXL,  L-BIP,  R-BIP  ) 
orient-«ith-PD(  L-FOOT-ORIEITATIOI ,  L-FOOT,  foot-stnfce-orientation,  tiae-m-air  ) 
set-angle-eith-PD!  L-TOI-AIGLE,  bL-TOI-JOIIT,  180®,  tiae-ia-air  ) 

Cases  L-SWING,  1-1ANDING,  and  R- TAKEOFF 
are  aa elogous  to  the  precedinq  three  cases. 

and 

cad 


Figure  9:  Abbreviated  Walking  Algorithm 
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